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Abstract 


HVACSIK*^  is  a modular,  non-proprietary  computer  simulation  package  de- 
veloped at  the  National  Bureau  of  Standards,  designed  to  allow  detailed  simu- 
lation of  entire  building  systems:  the  heating,  ventilating,  and  air  condi- 
tioning (HVAC)  system,  the  equipment  control  system,  the  building  shell,  the 
physical  plant,  and  the  dynamic  interactions  among  these  subsystems. 

The  HVACSIM'*'  package  consists  of  a main  simulation  program,  a library 
of  subroutines  containing  mathematical  models  of  building  energy  system  compo- 
nents, and  two  programs  used  in  preparing  a description  of  the  system  to  be 
simulated.  Models  representing  the  components  of  a physical  plant,  such  as 
boilers  and  chillers,  and  a model  representing  a multizone  building,  are  under 
development  and  will  be  added  to  the  HVACSIM'*’  package  as  they  become  avail- 
able. 


This  document  describes  the  operation,  capabilities,  and  limitations  of 
the  main  simulation  program,  presents  descriptions  of  the  existing  library  of 
component  models  for  simulating  HVAC  systems  and  their  controls,  includes 
theoretical  or  experimental  validation  of  the  models  wherever  possible,  and 
explains  the  procedure  by  which  new  component  models  may  be  written  and  added 
to  the  library. 


KEY  WORDS:  building  dynamics;  building  simulation;  building  system  modeling; 
computer  simulation  techniques;  control  dynamics;  dynamic  modeling  of  building 
systems;  dynamic  performance  of  building  systems;  dynamic  simulation;  HVAC 
system  simulation;  HVACSIM+ 
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INTRODUCTION  AND  OVERVIEW 


1 . 


1.0  INTRODUCTION 

HVACSIM'*’  is  a non-propr ietary  computer  simulation  package  developed  at 
the  National  Bureau  of  Standards  (NBS),  designed  to  allow  detailed  simulation 
of  entire  building  systems:  the  heating,  ventilating,  and  air  conditioning 
(HVAC)  system,  the  equipment  control  system,  the  building  shell,  the  physical 
plant,  and  the  dynamic  interactions  among  these  subsystems. 

The  HVACSIM"^  package  consists  of  a main  simulation  program  called  MODSIM, 
a library  of  subroutines  containing  mathematical  models  of  building  energy 
system  components,  an  interactive  'front  end'  program  called  HVACGEN  which  is 
used  to  create  a description  of  the  system  to  be  simulated,  and  a program 
called  SLIMCON  which  processes  the  output  of  HVACGEN  and  converts  it  into  the 
form  required  by  MODSIM. 

MODSIM  by  itself  is  a general-purpose  simulation  program,  employing  ad- 
vanced equation  solving  techniques  and  a hierarchical,  modular  approach  to 
solve  large  sets  of  non-linear  differential  and  algebraic  equations.  In  prin- 
ciple, MODSIM  could  be  used  to  simulate  any  system  which  can  be  represented  as 
a set  of  discrete,  interconnected  components.  The  inclusion  of  a library  of 
subroutines  describing  HVAC  system  components  makes  MODSIM  specific  to  build- 
ing energy  systems.  Models  representing  the  components  of  a physical  plant, 
such  as  boilers  and  chillers,  and  a model  representing  a multi-zone  building, 
are  under  development  and  will  be  added  to  the  HVACSIM^  package  as  they  become 
available. 

This  document  is  intended  to  serve  as  a reference  manual  for  HVACSIM^.  A 
separate  document  [1]  describes  the  operation  of  HVACGEN  and  the  installation 
and  use  of  the  HVACSIM^  package,  and  provides  a series  of  examples  to  serve  as 
a tutorial  in  the  use  of  HVACSIM^. 

The  remainder  of  this  chapter  describes  the  operation,  capabilities,  and 
limitations  of  MODSIM.  Most  of  this  information  has  been  drawn  from  reference 
[2],  which  was  written  by  the  author  of  the  MODSIM  program.  A brief  summary 
of  the  entire  simulation  package  is  also  provided. 

Chapter  two  explains  the  interface  between  MODSIM  and  the  individual  com- 
ponent models.  The  purpose  of  this  chapter  is  to  describe  how  component  mod- 
els may  be  written  and  added  to  the  library. 

Chapter  three  describes  a number  of  functions  and  subroutines  which  are 
used  by  some  of  the  component  models  and  are  available  for  use  by  future  com- 
ponent models.  These  routines  represent  transport  delays,  hysteresis  effects, 
Bessel  functions,  linear  least  squares  curve  fitting,  dry  fin  efficiencies, 
and  the  properties  of  fluids. 

Chapter  four  provides  descriptions  of  each  of  the  component  models  pre- 
sently included  in  HVACSIM^.  Wherever  possible,  theoretical  or  experimental 
verification  of  the  models  has  been  included. 


1 


1.1  AN  OVERVIEW  OF  MODSIM 


1.1.1  Types,  Units.  Blocks,  and  Snperblocks 

MODSIM  stands  for  MODular  SIMulation.  Many  ideas  for  the  design  of  MODSIM 
were  borrowed  from  the  transient  systems  simulation  program  TRNSYS,  which  was 
developed  at  the  University  of  Wisconsin  Solar  Energy  Laboratory  [3].  Both 
programs  use  separate  subroutines,  called  TYPE  subroutines,  to  represent  each 
kind  of  system  component.  When  a simulation  is  constructed,  the  components 
modeled  by  the  TYPE  subroutines  are  linked  together.  Two  or  more  components  of 
the  same  kind,  such  as  two  pumps,  are  modeled  by  the  same  TYPE  subroutine. 
The  distinction  between  components  of  the  same  kind  is  maintained  by  assigning 
a UNIT  number  to  each  component  in  the  simulation.  Each  unit  has  distinct  in- 
puts, outputs,  parameters,  and  saved  workspace.  As  far  as  a user  is  con- 
cerned, a unit  can  be  viewed  as  a black  box  which  takes  a set  of  inputs  and 
produces  a set  of  outputs. 

A major  difference  between  MODSIM  and  many  other  building  simulation  pro- 
grams is  the  method  used  to  obtain  a self-consistent  solution  at  each  time 
step.  Instead  of  the  commonly  used  method  of  subsequent  substitution,  MODSIM 
uses  a simultaneous  nonlinear  equation-solving  package  called  SNSQ  [4],  which 
is  based  on  a modification  of  the  Powell  hybrid  method  [5]. 

Large  scale  simulations  require  simultaneous  solution  of  a large  number 
of  algebraic  and  differential  equations.  When  the  number  of  equations  is  much 
greater  than  about  twenty  or  thirty,  direct  solution  of  the  equation  set  using 
a program  such  as  SNSQ  becomes  prohibitively  time-consuming.  MODSIM  uses  a 
hierarchical  approach  to  this  problem,  by  partitioning  large  sets  of  equations 
into  smaller  subsets.  Between  the  entire  simulation  and  the  individual  units, 
two  intermediate  levels  of  structure  are  defined:  BLOCKS  and  SUPERBLOCKS. 

A block  is  simply  a set  of  units,  more  or  less  arbitrarily  defined  by  the 
user.  Usually,  blocks  should  be  chosen  to  represent  functional  units  of 
moderate  size,  such  as  an  air  handler  or  the  elements  of  a control  loop.  The 
interconnections  among  the  units  in  a block  define  a set  of  simultaneous 
equations  which  is  solved  by  the  equation  solver  SNSQ.  Some  of  the  inputs  to 
the  units  in  a block  will  be  solved  simultaneously  within  the  block.  Any  unit 
inputs  which  are  not  solved  simultaneously  are  called  'block  inputs.'  These 
block  inputs  can  be  viewed  as  boundary  conditions  for  the  set  of  simultaneous 
equations  being  solved  within  the  block.  All  the  outputs  of  the  units  in  a 
block  can  also  be  viewed  as  block  outputs.  The  values  of  the  block  outputs 
will  depend  on  the  values  of  the  block  inputs.  We  can  now  view  a block  as  a 
black  box  which  takes  a set  of  block  inputs  and  produces  a set  of  block 
output  s . 

The  next  level  of  structure  is  the  superblock,  which  consists  of  a set  of 
blocks.  Blocks  within  a superblock  are  exactly  analogous  to  units  within  a 
block.  Connections  among  the  blocks  in  a superblock  define  a set  of  simulta- 
neous equations  which  is  solved  by  SNSQ.  A superblock  forms  a functional  unit 
of  a building  system,  and  a simulation  may  consist  of  several  superblocks. 
For  example,  one  superblock  might  contain  several  air  handlers,  another  might 
consist  of  air  conditioning  zones  within  a building,  and  a third  might  repre- 
sent the  building  shell.  However,  MODSIM  assumes  that  the  couplings  among  su- 
perblocks are  weak  enough  so  that  each  superblock  can  be  treated  as  an  inde- 
pendent subsystem.  No  simultaneous  equations  are  defined  between  superblocks. 
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1.1.2  Time  Steps  and  the  Integration  of  Differential  Equations 

Differential  equations  limit  the  size  of  time  steps  which  can  be  taken 
during  the  course  of  a simulation.  If  the  time  step  is  too  large,  inaccuracy 
or  instability  will  result.  This  is  particularly  a problem  when  the  compo- 
nents of  a system  have  a wide  range  of  characteristic  times,  or  time  con- 
stants, as  is  generally  the  case  in  building  energy  systems.  A set  of  differ- 
ential equations  with  widely  varying  time  constants  is  referred  to  as  a stiff 
set  of  equations. 

The  best  integration  algorithms  for  dealing  with  stiff  equations  are  the 
so-called  variable  time  step,  variable  order  Gear  algorithms  [6,7].  These  al- 
gorithms make  use  of  backward  differentiation  formulas  and  the  predictor-cor- 
rector method.  The  order  of  integration  and  the  size  of  the  time  step  are 
varied  during  the  calculation  to  minimize  the  computation  required  for  a spe- 
cified degree  of  accuracy.  The  algorithm  implemented  in  MODSIM  is  a modifica- 
tion of  the  algorithm  described  by  Brayton  [7]. 

Each  superblock  in  a simulation  is  an  independent  subsystem  in  the  sense 
that  it  proceeds  forward  in  time  independently.  The  time  step  and  integration 
order  are  determined  independently  for  each  superblock.  Thus,  for  example, 
rapid  control  actions  in  one  superblock  may  cause  that  superblock  to  require 
very  short  time  steps,  but  the  time  step  taken  by  other  superblocks  can  be 
much  longer. 

1.1.3  Time  Dependent  Boundary  Conditions 

In  the  context  of  MODSIM,  a boundary  condition  or  boundary  variable  is 
defined  as  a state  variable  which  is  external  to  the  system  being  simulated, 
in  the  sense  that  it  is  not  an  output  of  any  component  of  the  system.  A typi- 
cal example  is  the  outside  air  temperature.  Boundary  variables  may  be  con- 
stant or  time  dependent.  Any  number  of  boundary  variables  may  be  specified  as 
time  dependent,  and  their  values  are  stored  in  a data  file.  As  a simulation 
progresses,  data  are  read  from  this  file  and  used  to  change  the  appropriate 
state  variables. 

The  first  column  of  the  boundary  variable  data  file  must  contain  values 
of  time.  Time  intervals  need  not  be  equal.  MODSIM  uses  third  order  Lagrangi- 
an  interpolation  to  find  boundary  values  at  the  necessary  times.  The  remain- 
ing columns  contain  values  of  the  boundary  variables  at  the  specified  times. 

MODSIM  begins  all  simulations  at  time  zero.  If  the  first  time  in  the 
boundary  variable  file  is  3600  seconds,  for  example,  MODSIM  will  still  simu- 
late the  time  interval  from  0 to  3600  seconds.  During  this  time  interval,  the 
boundary  variable  values  specified  for  the  first  time  in  the  boundary  file 
will  be  used.  To  avoid  confusion,  it  is  recommended  that  boundary  data  files 
always  begin  with  time  zero. 

In  certain  cases,  the  change  in  a boundary  variable  may  be  discontinuous. 
A typical  example  is  a set  point  change.  In  such  a situation  it  is  necessary 
to  reset  the  differential  equation  integration  algorithm  at  the  time  of  the 
discontinuity.  This  is  accomplished  by  including  two  data  records  with  the 
same  time.  When  MODSIM  encounters  two  identical  times  in  sequence  in  the 
boundary  file,  it  integrates  up  to  the  repeated  time  and  then  reduces  the 
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simulation  time  step  to  the  minimum  and  reduces  the  integration  algorithm  to 
first  order,  beginning  at  the  repeated  time.  This  reset  convention  may  be 
used  to  mark  abrupt  changes  in  any  of  the  boundary  conditions.  There  is  no 
limit  on  the  number  of  resets  allowed  in  a simulation,  but  at  least  three 
records  must  occur  between  resets. 

1 .1 .4  Initialization 

Each  time  MODSIM  is  run,  a complete  description  of  the  final  state  of  the 
simulation  is  written  to  a file.  An  initialization  option  has  been  included 
which  allows  the  final  state  of  one  simulation  to  be  used  as  the  initial  state 
of  another.  This  feature  has  two  obvious  uses.  First,  a series  of  simula- 
tions of  a given  system  can  be  run  from  the  same  initial  steady  state,  without 
having  to  repeat  a computationally  expensive  startup  transient  for  each  run. 
Second,  a long  simulation  can  be  run  in  several  conveniently  short  pieces, 
with  each  piece  beginning  precisely  where  the  previous  piece  ended. 

As  discussed  in  the  previous  section,  MODSIM  begins  all  simulations  at 
time  zero.  Thus  if  a long  simulation  is  to  be  run  in  several  pieces,  it  will 
generally  be  necessary  to  prepare  a separate  boundary  variable  file  (beginning 
at  time  zero)  for  each  piece. 

1.1.5  Report  Generation 

Two  kinds  of  report  are  generated  by  MODSIM:  a 'data  file'  containing  raw 
data  at  each  time  step  taken  during  the  simulation,  and  a 'list  file'  contain- 
ing a summary  of  the  simulation  configuration,  any  diagnostics  generated 
during  the  computation,  and  an  interpolated  listing  of  the  reported  variables 
at  equal  time  intervals.  A reporting  interval  and  a set  of  reported  variables 
may  be  specified  for  each  superblock.  Reports  are  generated  for  each  super- 
block  individually  as  the  simulation  progresses.  As  discussed  in  the  previous 
section,  an  output  file  containing  the  final  state  of  the  simulation  is  also 
written. 

1.1.6  State  Variable  Freezing  and  Block  Inactivation 

During  the  course  of  a simulation,  some  of  the  state  variables  may  reach 
steady  state,  i.e.  cease  to  vary  with  time.  If  such  a state  variable  is  being 
solved  simultaneously,  a saving  of  computation  time  can  be  realized  by  remov- 
ing it  from  the  set  of  simultaneous  equations.  This  is  referred  to  as  'freez- 
ing' the  variable.  After  a variable  is  frozen,  it  must  be  monitored  so  that 
it  can  be  returned  to  the  calculation,  or  'unfrozen,'  if  necessary. 

If  all  the  simultaneous  equations  in  a block  are  frozen,  and  all  its 
block  inputs  are  frozen,  the  block  can  be  marked  inactive.  Once  a block  is 
inactivated,  it  is  no  longer  necessary  to  monitor  the  frozen  state  variables 
in  the  block.  A block  is  marked  active  as  soon  as  one  of  its  block  inputs 
becomes  unfrozen. 

1.2  LIMITATIONS  OF  MODSIM 

All  simulation  programs  have  inherent  limitations,  and  MODSIM  is  no 
exception.  In  section  1.1,  a number  of  techniques  used  to  increase  the  compu- 
tational efficiency  of  MODSIM  over  and  above  that  of  brute  force  calculation 
were  described.  Not  surprisingly,  these  techniques  involve  tradeoffs.  This 
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section  discusses  the  drawbacks  of  the  measures  taken  to  achieve  the  goal  of 
fast  execution  speed  for  very  large  scale  simulations,  and  describes  some 
additional  measures  taken  to  minimize  the  difficulties. 

1.2.1  Convergence  of  the  Equation  Solver 

Convergence  is  particularly  a problem  at  the  beginning  of  a simulation, 
when  the  equation  solver  must  work  from  the  set  of  user-supplied  initial 
conditions  chosen  somewhat  arbitrarily  to  start  the  calculation.  If  the 
initial  conditions  (or  an  intermediate  solution)  are  too  far  from  the  final 
solution,  the  equation  solver  may  fail  to  converge.  In  some  cases,  errors 
made  by  the  equation  solver  can  cause  erratic  solutions  for  several  time 
steps,  even  if  the  solutions  are  convergent.  The  choice  of  the  error  toler- 
ance for  the  equation  solver  is  an  important  factor  which  influences  this 
behavior.  For  further  discussion  of  error  tolerances,  see  section  1.3  below. 
The  convergence  properties  of  the  equation  solver  are  also  strongly  influenced 
by  the  block/superblock  structure  defined  when  a simulation  is  created. 

The  problem  of  slow  and  uncertain  convergence  to  a self-consistent  set  of 
initial  conditions  is  unavoidable,  but  the  consequences  can  be  minimized  by 
the  use  of  the  initialization  option.  This  permits  a series  of  simulations  to 
be  run  from  the  same  initial  state,  without  repeating  an  inaccurate  and  compu- 
tationally expensive  startup  transient  at  the  beginning  of  each  run. 

1.2.2  Limitations  of  Block  Structure 

It  was  stated  above  that  the  convergence  properties  of  the  equation 
solver  depend  on  the  block/superblock  structure.  This  is  not  surprising  since 
the  bl ock/ supe rbl ock  structure  determines  the  equation  sets  to  be  solved 
during  a simulation.  This  fact  limits  the  flexibility  of  simulation  construc- 
tion. Some  care  must  be  exercised  in  defining  blocks  to  obtain  good  conver- 
gence properties,  particularly  when  control  loops  are  involved.  In  general, 
it  is  good  practice  to  have  a control  loop  entirely  contained  in  a single 
block. 

1.2.3  Coupling  Between  Superblocks 

Superblocks  are  assumed  to  be  only  weakly  coupled.  No  simultaneous  equa- 
tions are  defined  between  them,  and  they  are  allowed  to  evolve  independently 
in  time.  Consequently,  great  care  is  required  in  dividing  a system  into 
superblocks. 

When  two  superblocks  are  coupled,  information  is  passed  between  them.  If 
the  coupling  between  superblocks  involves  a control  loop,  or  variables  which 
ideally  should  be  solved  simultaneously  such  as  pressures  and  flow  rates,  the 
calculations  tend  to  be  inaccurate  or  unstable.  To  minimize  this  problem, 
special  purpose  subroutines  have  been  designed  for  passing  information  (con- 
trol signals  and  flow  streams)  between  superblocks.  However,  the  weak  cou- 
pling of  superblocks  remains  a fundamental  limitation. 

Difficulties  can  also  arise  due  to  the  independent  time  evolution  of 
superblocks,  since  the  inputs  to  a superblock  may  change  at  a time  when  that 
superblock  is  not  scheduled  to  be  called.  To  mitigate  this  problem,  a super- 
block input  scanning  option  (INSOPT)  is  available.  When  this  option  is  se- 
lected (INSOPT  = 1),  all  superblock  inputs  are  scanned  after  each  time  step. 
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If  the  inputs  to  a superblock  which  was  not  called  during  the  time  step  have 
changed  by  an  amount  greater  than  an  error  tolerance,  the  superblock  is  called 
and  its  state  is  calculated  based  on  the  new  values  of  its  inputs.  This 
option  increases  both  the  accuracy  and  the  computational  requirements  of  a 
simulation.  When  the  input  scanning  option  is  not  selected  (INSOPT  = 0),  each 
superblock  is  called  at  times  determined  by  the  time  step  control  algorithm, 
regardless  of  any  changes  in  its  inputs  at  intermediate  times. 

1.2.4  Time  Step  Control 


Under  certain  conditions,  the  variable  time  step  integration  algorithm 
can  produce  erratic  results.  For  the  most  part,  these  problems  are  eliminated 
by  incorporating  a mechanism  for  rejecting  time  steps  which  are  too  large,  and 
by  defining  a minimum  and  a maximum  time  step.  The  choice  of  minimum  and 
maximum  time  steps  influences  the  accuracy  and  stability  of  the  simulation,  as 
well  as  the  amount  of  computation  required. 


The  time  step  selection  algorithm  depends  on  differential  equation  inte- 
gration errors  for  its  operation.  A simulation  which  contains  no  differential 
equations  will  always  proceed  at  the  maximum  time  step.  If  all  differential 
equations  in  a superblock  become  frozen  during  the  course  of  a simulation,  the 
length  of  each  time  step  taken  by  the  superblock  is  twice  the  length  of  the 
previous  time  step,  until  the  maximxim  step  size  is  reached. 


1.2.5  State  Variable  Freezing  and  Unfreezing 

A state  variable  is  frozen  if  it  changes  less  than  a certain  amount, 
determined  by  user-supplied  error  tolerances,  from  one  time  step  to  the  next. 
Clearly,  state  variable  freezing  will  introduce  some  errors  into  the  computa- 
tion. When  a state  variable  is  unfrozen  within  a block,  the  block  equations 
are  solved  again  with  the  unfrozen  variable  put  back  into  the  equation  set. 
If  this  causes  any  additional  variables  to  unfreeze,  all  of  the  variables 
solved  simultaneously  within  the  block  are  unfrozen  and  the  computation  is 
repeated  a third  time.  Thus  the  unfreezing  of  variables  forces  extra  computa- 
tion, both  by  increasing  the  size  of  the  equation  set  and  by  causing  calcu- 
lations to  be  repeated. 


For  superblock  equations,  three  optional  modes  are  implemented.  The  mode 
is  determined  by  the  value  of  the  'superblock  variable  unfreezing  option,' 
IFZOPT.  If  a superblock  equation  is  unfrozen  and  mode  0 is  set,  the  corre- 
sponding state  variable  is  not  put  back  into  the  equation  set  until  the  next 
time  step.  If  mode  1 is  set,  any  unfrozen  variables  are  put  back  into  the 
superblock  equation  set  and  the  calculation  is  repeated.  If  mode  2 is  set, 
all  superblock  equations  are  put  back  into  the  equation  set  and  the  calcula- 
tion is  done  again.  The  reason  for  these  three  modes  is  that  superblock 
equations  require  much  more  computation  than  block  equations.  In  some  cases, 
depending  on  the  block/superblock  structure,  the  errors  introduced  by  choosing 
mode  0 or  1 are  acceptable,  and  these  modes  result  in  a significant  saving  of 
computation  time.  In  other  cases  the  simpler  modes,  particularly  mode  0, 
result  in  serious  errors,  and  mode  3 is  required  for  satisfactory  accuracy. 

1.2.6  Size  of  Simulations 


The  size  of  simulations  is 
ty  of  arrays  in  MODSIM.  These 


limited  by  the  dimensions  assigned  to  a varie- 
are  not  fundamental  limitations,  since  dimen- 
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sion  statements  can  be  changed.  As  distributed,  MODSIM  has  the  following 
1 imits : 

15  inputs  per  unit 

15  outputs  per  unit 

10  differential  equations  per  unit 

20  units  per  block 

50  inputs  per  block 

50  outputs  per  block 

30  simultaneous  equations  between  units  in  a block 
10  blocks  per  superblock 

20  differential  equations  per  superblock 

20  simultaneous  equations  between  blocks  in  a superblock 
30  reported  variables  per  superblock 
200  units  per  simulation 
50  blocks  per  simulation 
10  superblocks  per  simulation 
50  differential  equations  per  simulation 
600  state  variables  per  simulation 
1000  parameters  per  simulation 
1000  SAVED  variables  per  simulation 

30  time-dependent  boundary  variables  per  simulation 

1,3  AN  OVERVIEW  OF  HVACSIM"^ 

This  section  provides  a brief  summary  of  the  procedure  for  preparing  and 
running  simulations  using  HVACSIM^.  A more  complete  and  detailed  account  can 
be  found  in  reference  [1]. 

First  it  is  necessary  to  decide  how  the  system  to  be  simulated  will  be 
divided  into  units,  blocks,  and  superblocks.  A diagram  of  the  system  is 
indispensable  at  this  stage.  TYPE  subroutines  representing  each  unit  must  be 
selected,  all  state  variables  must  be  identified  and  numbered,  and  values  for 
the  parameters  required  by  each  component  must  be  determined. 

A description  of  the  system  to  be  simulated  is  then  prepared  using  the 
interactive  program  HVACGEN.  Information  pertaining  to  the  configuration  of 
each  TYPE  component  is  contained  in  a file  called  TYPAR.DAT,  which  is  read  by 
HVACGEN.  HVACGEN  prompts  the  user  to  specify  the  TYPE  of  each  component  in 
the  simulation,  to  provide  index  numbers  for  each  state  variable,  and  to  enter 
values  for  all  parameters.  HVACGEN  assigns  a UNIT  number  to  each  component, 
to  distinguish  between  components  of  the  same  type. 

After  all  UNITS  have  been  entered,  HVACGEN  asks  the  user  to  supply  sever- 
al error  tolerances: 

RTOLX,  ATOLX  - relative  and  absolute  error  tolerances.  These  tolerances  are 
used  in  integrating  differential  equations  and  choosing  the  time  step 
size,  and  in  deciding  when  to  freeze  or  unfreeze  state  variables.  Some 
TYPE  subroutines  and  utility  functions  also  use  these  tolerances  to 
decide  when  internal  iterative  procedures  have  converged  sufficiently, 
and  to  decide  when  the  freezing  of  differential  equations  should  be  pei^ 
mitted. 
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XTOL  - an  error  tolerance  used  by  the  equation  solver,  SNSQ.  The  value  of 
XTOL  influences  the  convergence  properties  of  the  equation  solver  (see 
section  1.2.1).  As  a rule  of  thumb,  XTOL  should  be  greater  than  or  equal 
to  the  sum  of  RTOLX  and  ATOLX. 

TTIME  - the  time  interval  over  which  the  integral  of  a differential  equation 
must  satisfy  the  error  criterion  defined  by  RTOLX  and  ATOLX. 

HVACGEN  supplies  somewhat  arbitrary  default  values  for  these  error  tolerances. 
The  default  tolerances  are  probably  unnecessarily  stringent  for  most  purposes, 
and  larger  values  could  be  used.  However,  it  is  important  that  XTOL  be 
greater  than  or  equal  to  the  sum  of  RTOLX  and  ATOLX.  In  other  words,  the 
criterion  for  freezing  and  unfreezing  variables  must  be  more  stringent  than 
the  convergence  criterion  for  the  equation  solver.  Otherwise  variables  may  be 
frozen  and  unfrozen  frequently  and  erratically,  resulting  in  a slow  and  unsta- 
ble simulation. 

Next  HVACGEN  asks  for  initial  values  of  all  state  variables,  indices  of 
time  dependent  boundary  variables,  a reporting  interval  for  each  superblock, 
and  indices  of  variables  to  be  reported. 

It  then  prompts  the  user  to  select  a superblock  variable  unfreezing  op- 
tion (IFZOPT)  and  a superblock  input  scan  option  (INSOPT).  The  variable 
freezing  option  determines  the  manner  in  which  superblock  equations  are  recal- 
culated when  a superblock  equation  unfreezes.  Three  modes  are  available,  as 
described  in  section  1.2.5.  The  value  of  IFZOPT  influences  the  strength  of 
coupling  among  the  blocks  within  a superblock,  and  has  no  effect  on  simula- 
tions consisting  of  a single  block.  The  input  scan  option  is  mentioned  in 
section  1.2.3.  Enabling  this  option  (mode  1)  allows  a superblock  to  be  called 
whenever  its  inputs  change,  overriding  the  time  step  selection  algorithm. 
When  the  input  scan  option  is  disabled  (mode  0),  a superblock  is  called  only 
at  times  determined  by  the  time  step  selection  algorithm,  regardless  of  any 
changes  in  the  superblock  inputs  at  intermediate  times.  The  value  of  INSOPT 
influences  the  coupling  between  superblocks,  and  has  no  effect  on  simulations 
consisting  of  a single  superblock. 

HVACGEN  also  has  extensive  editing  capabilities  to  facilitate  the  cre- 
ation, testing,  and  modification  of  large  system  simulations.  The  operation 
and  features  of  HVACGEN  are  the  subject  of  a separate  report  [1].  The  output 
of  HVACGEN  is  a 'simulation  work  file.' 

Before  a simulation  can  be  run,  the  program  SLIMCON  must  be  run.  This 
program  reads  the  simulation  work  file,  performs  some  processing  steps  on  the 
numbers,  and  writes  a 'model  definition  file'  which  contains  a complete  system 
description  in  the  form  required  by  MODSIM. 

When  the  simulation  program  is  run,  it  prompts  the  user  for  the  following 
information: 

IPRINT  - a non-zero  value  of  IPRINT  causes  intermediate  results  of  the  equa- 
tion solver  to  be  written  to  the  output  'list  file.'  This  feature  can  be 
useful  for  debugging,  but  produces  vast  quantities  of  output.  Normally 
IPRINT  is  set  to  zero,  and  only  final  values  at  each  reporting  time  are 
written. 
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INIT  ~ a non-zero  value  of  INIT  causes  the  initial  state  of  the  system  to  be 
read  from  a file.  If  INIT=0,  the  initial  state  vector  in  the  model 
definition  file  defines  the  initial  state  of  the  system. 

In  addition,  the  program  asks  for  a minimum  time  step,  a maximum  time 
step,  and  a stopping  time.  The  choice  of  minimum  and  maximum  time  steps  in- 
fluences the  accuracy  and  stability  of  the  simulation,  as  well  as  the  amount 
of  computation  required. 

As  the  simulation  progresses,  MODSIM  calls  the  TYPE  subroutines  repre- 
senting the  components  of  the  system.  Some  of  the  TYPE  subroutines  call  util- 
ity functions  and  subroutines  described  in  chapter  3.  The  architecture  of 
HVACSIM'*^  is  summarized  in  figure  1.3.1. 


Figure  1.3.1  Architecture  of  the  HVACSIM''"  package 
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2 .0  FORMULATION  OF  COMPONENT  MODELS 


The  library  of  component  models  described  in  this  report  includes  most  of 
the  components  of  HVAC  systems,  including  controllers  for  such  systems.  How- 
ever, the  library  is  by  no  means  exhaustive,  and  situations  will  inevitably 
arise  in  which  additional  component  models  are  needed. 

In  addition,  it  is  worth  mentioning  that  nothing  in  the  main  simulation 
program,  MODSIM,  limits  it  to  building  energy  analysis.  MODSIM  is  a general- 
purpose  modular  simulation  program  which  can  be  used  to  study  any  dynamic 
system  for  which  suitable  component  models  are  available.  It  is  assumed  that 
MODSIM  will  be  used  primarily  as  part  of  HVACSIM^,  by  researchers  involved  in 
the  analysis  of  HVAC  systems  and  controls,  though  other  applications  of  MODSIM 
are  possible. 

This  section  presents  general  information  on  the  structure  of  component 
model  TYPE  subroutines  consistent  with  the  requirements  of  MODSIM.  Develop- 
ment of  mathematical  models  to  represent  a physical  process  is  often  a diffi- 
cult undertaking,  requiring  a certain  degree  of  expertise.  However,  this 
problem  is  not  addressed  here.  Rather,  emphasis  is  on  features  influencing 
the  interface  between  TYPE  subroutines  and  MODSIM.  The  TYPE  9 linear  valve 
model  described  in  section  4.9  is  used  as  an  example  to  point  out  important 
features  of  TYPE  subroutines  in  general.  A listing  of  the  FORTRAN  source  code 
for  the  TYPE  9 subroutine  appears  at  the  end  of  section  2.2,  and  is  referred 
to  throughout  this  section.  The  reader  should  review  section  4.9  before 
proceeding  through  the  example.  Familiarity  with  FORTRAN  programming  is 
assumed.  HVACSIM^  (including  the  TYPE  subroutines)  is  written  in  ANSI  Stan- 
dard X3. 9-1978  FORTRAN,  which  is  compatible  with  FORTRAN  VII  compilers. 

2.1  SOLUTION  OF  DIFFERENTIAL  EQUATIONS 

The  model  definition  file  informs  MODSIM  of  the  number  of  differential 
equations  in  each  UNIT  of  the  simulation.  Provision  has  been  made  in  the 
setup  program,  HVACGEN,  to  allow  the  number  of  differential  equations  in  a 
given  component  to  be  determined  by  the  values  of  one  or  more  parameters. 
Thus  different  UNITS  of  a given  TYPE  may  have  different  numbers  of  differen- 
tial equations.  This  feature  is  discussed  further  in  section  2.3. 

For  a UNIT  with  N differential  equations,  MODSIM  expects  the  first  N 
outputs  to  be  first  derivatives  of  the  N dependent  variables,  rather  than  the 
dependent  variables  themselves.  The  dependent  variables  may  be  inputs  to  the 
component,  and  to  any  other  component  in  the  system.  No  limiting  assumptions 
are  made  about  the  form  of  differential  equations.  The  outputs  must  be  first 
order  derivatives,  but  this  is  not  a fundamental  limitation  since  higher  order 
differential  equations  can  be  factored  into  a set  of  coupled  first  order 
differential  equations. 

2.2  EXAMPLE:  LINE-BY-LINE  DESCRIPTION 
Line  1.: 

1 SUBROUTINE  TYPE9 (XIN, OUT, PAR, SAVED, lOSTAT) 

The  symbolic  name  of  each  TYPE  subroutine  is  of  the  form  TYPEn,  where  n 
is  an  integer  less  than  or  equal  to  99.  TYPE  subroutines  are  called  by 
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subroutine  SELECT,  which  must  include  a CALL  statement  for  each  TYPE  subrou- 
tine. The  arguments  of  the  subroutine  are  defined  as  follows: 


XIN  : 

OUT  : 

PAR 

SAVED  : 
lOSTAT  : 


an  array  of  time-dependent  inputs  to  the  component 

an  array  of  time-dependent  outputs  from  the  component 

an  array  of  constant  parameters  defining  the  characterisics  of  the 

component 

a saved  workspace  specific  to  a particular  implementation  or  UNIT  of 
a component 

an  array  of  information  related  to  the  freezing  and  unfreezing  of 
state  variables 


Lines  5-6 : 


5 REAL  K,XIN(4) ,0DT(2) ,PAR(4) ,SAVED(5) 

6 INTEGER  I0STAT(4) 

Since  the  arguments  of  the  subroutine  are  all  arrays,  they  must  be  dimen- 
sioned within  the  subroutine.  lOSTAT  has  the  same  dimension  as  either  XIN  or 
OUT,  whichever  is  larger.  The  use  of  lOSTAT  is  discussed  below.  SAVED  need 
not  be  dimensioned  if  the  saved  workspace  is  not  used. 

Lines  8-10 ; 


8 COMMON/ CHRONO/  TIME,TSTEP,TTIME,TMIN, ITIME 

9 COMMON/ SOS  COM/  RTOLX , ATOLX , XTOL 

10  COMMON/XINIT/  INIT,NSAVED 

MODSIM  uses  common  blocks  extensively  to  pass  information  internally. 

Three  common  blocks  are  invoked  by  the  TYPE9  subroutine. 

CHR(X40  contains  the  following  information: 

TIME  : the  current  simulation  time 

TSTEP  : the  current  simulation  time  step 

TTIME  : a time  interval  associated  with  error  criteria  for  integrating 
differential  equations 

TMIN  : the  minimum  allowable  simulation  time  step 

ITIME  : the  number  of  time  steps  taken  since  the  beginning  of  the  simula- 
tion 

SOSCOM  contains  the  following  information: 

RTOLX  : a relative  error  tolerance 

ATOLX  : an  absolute  error  tolerance,  used  together  with  RTOLX  to  calcu- 
late error  criteria  for  differential  equations  and  bounds  for 
variable  freezing 

XTOL  : an  error  tolerance  for  simultaneous  equation  solving 

XINIT  contains  the  following  information: 

INIT  : a flag  determining  whether  the  simulation  is  to  be  initialized 

from  the  initial  state  vector  (INIT=0)  or  from  an  initialization 
file  (INITIO) 

NSAVED:  the  number  of  saved  workspace  variables  in  the  entire  simulation 
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Lines  12-19: 


12  P2=XIN(1) 

13  W=XIN(2) 

14  C=XIN(3) 

15  CV=XIN(4) 

16  K=PAR(1) 

17  TC0N=PAR(2) 

18  CL=PAR(3) 

19  HYS=PAR(4) 

P2  is  the  water  pressure  at  the  valve  outlet,  W is  the  water  mass  flow 
rate,  C is  a control  signal  specifying  a requested  valve  position,  and  CV  is 
the  actual  position  of  the  actuator.  The  value  of  CV  is  obtained  from  the 
first  output  of  this  subroutine,  as  discussed  below  after  lines  25  and  39. 
The  four  parameters  are  the  flow  resistance  when  the  valve  is  open,  the  time 
constant  of  the  actuator,  a leakage  parameter,  and  a hysteresis  parameter, 
respectively.  Note  that  mnemonics  for  the  parameters  must  be  assigned  on  each 
call  to  the  subroutine,  even  though  they  do  not  change  with  time,  since  this 
subroutine  may  be  used  to  represent  two  or  more  valves  with  different  parame- 
ters in  the  same  simulation. 

Lines  21-24 : 


21  IF(C.LT.O.)  C=0. 

22  IF(C.GT.l.O)  C=1.0 

23  IF(CV.LT.O.)  CV=0. 

24  IF(CV.GT.l.O)  CV=1.0 

In  the  process  of  obtaining  a solution,  the  MODSIM  equation  solver  may 
occasionally  call  the  subroutine  with  values  of  C or  CV  greater  than  one  or 
less  than  zero.  These  lines  enforce  the  requirement  that  the  control  signals 
C and  CV  must  lie  between  0 and  1,  inclusive. 

Lines  25  , 39,  and  46 : 

25  IF(TCON.GE.l.O)  GOTO  1000 

39  1000  DCVDT=(C-CV) /TCON 

46  0UT(1)=DCVDT 

The  description  of  the  TYPE9  component  in  TYPAR.DAT  informs  MODSIM  that 
the  component  has  one  differential  equation  if  the  value  of  the  second  parame- 
ter is  greater  than  or  equal  to  one,  and  no  differential  equations  otherwise. 
Thus  if  the  time  constant  of  the  actuator  is  one  second  or  more,  lines  26-38 
are  skipped  and  the  dynamics  of  the  actuator  are  represented  by  line  39.  Since 
MODSIM  expects  one  differential  equation,  it  assumes  the  first  output  is  the 
time  derivative  of  a state  variable,  rather  than  the  state  variable  itself. 
MODSIM  takes  care  of  integrating  the  derivative  to  evaluate  the  state  vari- 
able, in  this  case  CV.  If  the  actuator  time  constant  is  less  than  one  second, 
the  differential  equation  is  solved  within  the  subroutine  and  the  first  output 
is  the  value  of  CV,  rather  than  its  derivative. 
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Lines  26-29  and  45 : 


26 

27 

28 

29  3000 

45  SAVED(4)=CV 

The  SAVED  array  is  generally  used  to  store  previous  values  of  variables 
whose  history  is  required.  Normally  two  spaces  in  the  SAVED  array  are  needed 
for  each  variable  to  be  saved,  plus  one  space  for  storing  the  simulation  time. 

In  this  particular  case,  SAVED(4)  holds  the  value  of  CV  from  the  previous 
call  to  the  UNIT,  which  may  be  called  repeatedly  at  a single  simulation  time. 
SAVED(5)  holds  the  final  value  of  CV  from  the  last  call  of  the  previous 
time  step,  while  SAVED(3)  holds  the  simulation  time  from  the  previous  call. 
Thus  line  45  saves  CV  on  each  call,  and  line  29  updates  SAVED(5)  only  on  the 
first  call  of  a new  timestep.  The  utility  function  HYSTER,  called  in  line  42, 
takes  care  of  maintaining  the  simulation  time  in  SAVED(3). 

Line  26  serves  to  identify  the  first  time  step  of  the  simulation.  If  an 
initialization  file  is  used  (INIT  0),  the  saved  workspace  initially  contains 
final  values  from  the  initialization  run.  In  this  case  SAVED(4)  should  not  be 
changed,  but  SAVED(3)  will  initially  contain  the  final  time  of  the  initializa- 
tion run,  and  must  be  reset  to  zero  on  the  first  call  of  the  first  time  step, 
as  is  done  in  line  28.  If  no  initialization  file  is  used,  both  SAVED(3)  and 
SAVED(4)  are  initialized  in  lines  27  and  28. 

Lines  30-38  and  46 : 


IF(ITIME.GT.l)  GOTO  3000 
IF ( INIT. EQ. 0) SAVED( 4) =CV 

IFdNIT.EQ.O  .OR.  SAVED (3 ). GT. TIME)  SAVED(3)=0. 
IF(TIME.GT.SAVED(3) ) SAVED( 5) =SAVED( 4) 


30  DCVDT=C 

31  CV=C 

32  A=TCON/TSTEP 

33  IF(A.LT.O.OS)  GOTO  2000 

34  DC=C-SAVED(5) 

35  IF(ABS(DC) .LT.l.E-10)  GOTO  2000 

36  CV=C-DC*EXP(-1./A) 

37  DCVDl^CV 

38  GOTO  2000 

46  OUT(l)=DCVDT 


If  the  actuator  time  constant  is  less  than  one  second,  the  differential 
equation  in  line  39  is  solved  in  lines  30-38.  The  resulting  value  of  CV  is 
assigned  to  the  local  variable  DCVDT  as  a matter  of  convenience,  allowing  line 
46  to  be  used  regardless  of  whether  the  differential  equation  is  solved  inside 
or  outside  the  subroutine.  Note  that  care  has  been  taken  to  avoid  division  by 
zero  when  the  time  constant  is  zero,  and  to  avoid  causing  FORTRAN  overflow  or 
underflow  errors  when  CV  is  very  nearly  equal  to  C. 
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Lines  40-43  and  46-41 : 


40  2000 

S=l. 

41 

IF(W.NE.0.)  S=SIGN(1.,W) 

42 

CH=HYSTER ( CV, SAVED ( 1 ) , SAVED ( 2 ) , SAVED ( 3 ) , HYS ) 

43 

P1=P2+S*K*(W/ ( (l.-CL)*CH+CL) )**2 

46 

0DT(1)=DCVDT 

47 

0UT(2)=P1 

These 

lines  complete  the  mathematical  description  of 

a linear  val 

hysteresis 

, as  described  in  Section  4.9,  and  place  the 

results  in  i 

array.  The 

utility  function  HYSTER  is  described  in  Chapter  3. 

Lines  7 and 

48-51: 

7 

LOGICAL  CON 

48 

CON=(IOSTAT(3) .LT.-l) .OR. (I0STAT(3) .EQ.2) 

49 

IOSTAT(1)=0 

50 

IF( (ABS(C-CV) .LE.RT0LX*ABS(CV)+AT0LX) .AND. CON) 

I0STAT(1)=1 

51 

IOSTAT(2)=l 

OUT 


lOSTAT  is  an  input/output  status  vector.  On  entry  this  vector  contains 
the  status  of  the  input  variables.  Possible  values  of  lOSTAT  and  their  mean- 
ings are  as  follows: 


lOSTAT  = 


-4  time  independent  boundary  condition 
-3  time  dependent  boundary  condition 
-2  frozen  output  (not  solved  simultaneously) 
-1  output  (not  solved  simultaneously) 

0 solved  simultaneously,  may  not  be  frozen 

1 solved  simultaneously,  may  be  frozen 

2 frozen 

3 unfrozen 


On  exit  it  contains  flags  which  enable  (I0STAT=1)  or  disable  (IOSTAT=0)  varia- 
ble freezing  for  the  output  variables.  In  general  it  is  recommended  that 
freezing  be  enabled  for  algebraic  outputs,  as  in  line  51.  However,  care  must 
be  taken  in  establishing  criteria  for  enabling  the  freezing  of  dependent 
variables  in  differential  equations.  In  line  48,  the  logical  variable  CON  is 
true  if  the  input  control  signal  C is  either  a frozen  variable  or  a boundary 
condition.  In  line  50,  freezing  of  the  dependent  variable  CV  is  enabled  only 
if  CON  is  true  and  CV  is  within  error  tolerances  of  C.  This  is  to  prevent  the 
differential  equation  from  freezing  if  CV  happens  to  be  nearly  equal  to  C at  a 
time  when  C is  changing  rapidly.  lOSTAT  should  always  be  set  to  either  0 or  1 
for  each  output. 
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1 

2 C 

3 C 

4 C 

5 

6 

7 

8 
9 

10 

11  C 

12 

13 

14 

15 

16 

17 

18 

19 

20  C 

21 
22 

23 

24 

25 

26 

27 

28 

29  3000 

30 

31 

32 

33 

34 

35 

36 

37 

38 

39  1000 

40  2000 

41 

42 

43 

44  C 

45 

46 

47 

48 

49 

50 

51 

52 

53 


SUBROUTINE  TYPE9(XIN, OUT, PAR, SAVED, lOSTAT) 

Linear  valve  model. 

REAL  K,XIN(4) ,OUT(2) ,PAR(4) ,SAVED(5) 

INTEGER  IOSTAT(4) 

LOGICAL  CON 

COMMON/ CHRONO/  TIME, TSTEP, TTIME, TMIN , ITIME 
COMMON/ SOS COM/  RTOLX , ATOLX , XTOL 
COMMON/XINIT/  INIT,NSAVED 

P2=XIN(1) 

W=XIN(2) 

C=XIN(3) 

CV=XIN(4) 

K=PAR(1) 

TCON=PAR(2) 

CL=PAR(3) 

HYS=PAR(4) 

IF(C.LT.O.)  C=0. 

IF(C.GT.l.O)  C=1.0 
IF(CV.LT.O.)  CV=0. 

IF(CV.GT.l.O)  CV=1.0 
IF(TCON.GE.l.O)  GOTO  1000 
IF ( ITIME. GT.l)  GOTO  3000 
IF ( INIT. EQ . 0 ) SAVED ( 4 ) =CV 

IFdNIT.EQ.O  .OR.  SAVED(3)  .GT.TIME)  SAVED(3)=0. 
IF(TIME.GT.SAVED(3))  SAVED( 5) =SAVED(4) 

DCVDI^C 

CV=C 

A=TCON/ TSTEP 
IF(A.LT.0.05)  GOTO  2000 
DC=C-SAVED(5) 

IF(ABS(DC) .LT.l.E-10)  GOTO  2000 
CV=C-DC*EXP(-1./A) 

DCVDT=CV 
GOTO  2000 
DCVDT=(C-CV) /TCON 
S=l. 

IF(W.NE.O.)  S=SIGN(1.,W) 

CH=HYSTER ( C V, SAVED ( 1 ) , SAVED ( 2 ) , SAVED ( 3 ) , HYS ) 

P1=P2+S*K*(W/ ((l.-CL)*CH+CL) )**2 

SAVED(4)=CV 

OUT(l)=DCVDT 

0UT(2)=P1 

CON=(IOSTAT(3) .LT.-l) .OR. (IOSTAT(3)  .EQ.2) 

IOSTAT(1)=0 

IF((ABS(C-CV) .LE.RTOLX*ABS(CV)+ATOLX) .AND. CON)  IOSTAT(l)=l 

IOSTAT(2)=l 

RETURN 

END 
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2.3  INSTALLATION  OF  COMPONENT  MODELS 


Before  a new  TYPE  subroutine  can  be  used,  a CALL  statement  for  the 
subroutine  must  be  added  to  subroutine  SELECT,  which  is  part  of  MODSIM.  In 
addition,  a description  of  the  component  must  be  added  to  the  file  TYPAR.DAT, 
which  is  read  by  HVACGEN.  This  section  describes  the  format  of  the  informa- 
tion required  by  HVACGEN.  The  TYPAR.DAT  entry  for  the  TYPE  9 linear  valve 
model  appears  at  the  end  of  this  section. 

The  first  line  must  have  an  asterisk  in  the  first  column.  This  signals 
HVACGEN  that  a new  TYPE  description  is  about  to  begin.  The  remainder  of  the 
line  is  ignored. 

The  second  line  contains  the  TYPE  number,  followed  by  a description  of 
the  component  in  quotation  marks. 

The  third  line  contains  five  integers,  specifying  the  number  of  SAVED 
variables,  differential  equations,  inputs,  outputs,  and  parameters,  respec- 
tively. In  the  case  of  the  linear  valve  model,  there  are  five  saved  vari- 
ables, one  differential  equation,  four  inputs,  two  outputs,  and  four  parame- 
ters . 


Next  comes  a set  of  lines  describing  the  outputs,  followed  by  a set  of 
lines  describing  the  inputs.  For  the  linear  valve  there  are  two  lines  for  the 
outputs  and  four  lines  for  the  inputs.  In  general  there  is  one  line  for  each 
output  and  one  line  for  each  input.  Each  of  these  lines  contains  an  integer 
followed  by  two  character  fields  enclosed  in  quotation  marks.  The  integer 
identifies  the  category  of  the  input  or  output.  Eight  categories  are  avail- 
able, as  follows: 

1 - Pressure 

2 - Flow  rate 

3 - Temperature 

4 - Control  signal  (or  other  fraction) 

5 - Rotation  rate 

6 - Energy 

7 - Power 

8 - Humidity 

The  first  character  field  provides  a brief  description  of  the  input  or  output. 
The  second  character  field  gives  recommended  or  required  units  for  the  vari- 
able. 

After  the  output  and  input  description  lines,  one  or  more  optional  lines 
may  be  used  to  modify  the  number  of  differential  equations.  Each  such  line 
consists  of  two  integers,  a relational  operator  in  quotation  marks,  and  a real 
number.  The  first  integer  is  a parameter  number.  The  second  integer  is  a 
number  to  be  subtracted  from  the  number  of  differential  equations  if  the  spe- 
cified condition  is  true.  The  truth  or  falsity  of  the  condition  is  determined 
by  applying  the  relational  operator  between  the  parameter  and  the  real  number. 
In  the  case  of  the  linear  valve  model,  the  line 

2 1 'LT'  1.0 

is  equivalent  to  the  FORTRAN  statement 
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IF  (PAR(2)  .LT.  1.0)  NDE  = NDE  - 1 


where  NDE  is  the  number  of  differential  equations. 

Any  of  the  six  FORTRAN  relational  operators  may  be  used:  'EQ',  'NE', 
'GT',  'GE',  'LT',  or  'LE'.  In  addition,  two  or  more  conditions  may  be  joined 
with  a logical  AND,  using  as  relational  operators  'EQ.AND',  'NE.AND',  and  so 
on.  The  second  integer  on  any  line  using  an  AND  operator  is  ignored.  For 
example,  the  two  lines 

3 0 'GT.AND'  0. 

3 2 'LE'  1. 

are  equivalent  to  the  FORTRAN  statement 

IF  ((PAR(3)  .GT.  0.) .AND.(PAR(3)  .LE.  1.))  NDE  = NDE  - 2 

Lines  modifying  the  number  of  differential  equations  are  optional  and  may 
be  omitted  entirely.  There  is  no  limit  to  the  number  of  lines  permitted,  or 
to  the  number  of  conditions  which  may  be  joined  with  logical  ANDs.  A number 
sign  (^)  in  column  one  tells  HVACGEN  not  to  expect  any  more  conditions.  The 
number  sign  is  required  whether  or  not  z.rsy  conditions  are  present. 

Finally,  one  line  for  each  parameter  is  required.  These  lines  consist  of 
the  parameter  number  and  a single  character  field  describing  the  parameter. 
The  set  of  parameter  description  lines  may  be  followed  by  any  number  of  blank 
lines  or  comment  lines,  which  will  be  ignored  by  HVACGEN. 


Listing  of  TYPAR.DAT  entry  for  the  TYPE  9 linear  valve  model 


*««**********«***««***  *«*»**  ******************  4>*****«*  **«***••«*** 

9 'LINEAR  VALVE  WITH  PNEUMATIC  ACTUATOR' 

5 14  2 4 

4 'ACTUATOR  RELATIVE  POSITION'  'DIMENSIONLESS' 

1  'INLET  WATER  PRESSURE'  'KPA' 

1 'OUTLET  WATER  PRESSURE'  'KPA' 

2 'WATER  MASS  FLOW  RATE'  'KG/ S' 

4 'CONTROL  VARIABLE  FROM  CONTROLLER'  'DIMENSIONLESS' 

4 'ACTUATOR  RELATIVE  POSITION  (SAME  AS  FIRST  OUTPUT)'  'DIMENSIONLESS' 
2 1 'LT'  1.0 

1 'FLOW  RESISTANCE  [0 .001/ (KG-M) ] ' 

2 'ACTUATOR  TIME  CONSTANT  [SEC] ' 

3 'LEAKAGE  PARAMETER  [DIMENSIONLESS] ' 

4 'HYSTERESIS  PARAMETER  [DIMENSIONLESS] ' 
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3.  UTILITY  FUNCTIONS  AND  SUBROUTINES 


3.0  Introduc  t ion 


This  chapter  describes  a number  of  functions  and  subroutines  which  are 
available  for  use  by  TYPE  components.  Sections  3.1  and  3.2  describe  two 
functions,  DELAY  and  HYSTER,  which  are  used  by  a number  of  component  models  to 
represent  transport  delays  and  hysteresis  effects,  respectively.  Section  3.3 
describes  four  subroutines  which  are  used  to  determine  heat  exchanger  fin 
efficiencies.  Three  of  these  subroutines,  which  compute  modified  Bessel 
functions  and  perform  polynomial  curve  fitting,  may  also  be  useful  for  other 
applications.  The  remaining  sections  summarize  functions  and  subroutines 
which  calculate  fluid  properties.  Most  of  these  are  not  used  by  any  of  the 
existing  component  models,  but  are  available  for  future  use. 

3.1  Transport  Delays 

General  Description 

Transport  delays  in  conduits  are  modeled  by  the  function  DELAY,  which 
takes  the  following  arguments: 


TIN  : 
TAU  : 
A : 

TOLD  : 
DXOLD  : 

TIMOLD: 


the  current  value  of  the  temperature  (or  other  quantity)  to  be  delayed 
the  transport  time,  based  on  the  current  flow  rate 

an  array  of  six  coefficients  for  a polynomial  representing  the  temper- 
ature distribution  in  the  conduit 
the  previous  value  of  TIN  (one  time  step  ago) 

the  previous  value  of  DX,  the  fractional  distance  the  fluid  moved  in  a 
time  step 

the  previous  simulation  time 


The  first  two  arguments,  TIN  and  TAU,  are  supplied  by  the  calling  subroutine. 
The  remaining  arguments  should  be  stored  in  the  saved  workspace  by  the  calling 
subroutine.  Thus  the  function  might  be  used  as  in  the  following  line: 

T = DELAY(T,TAU,SAVED(1) ,SAVED(7) ,SAVED(8) ,SAVED(9)) 


where  SAVED(l)  through  SAVED(6)  are  used  to  store  the  array  of  'A'  coeffi- 
cient s . 


Mathematical  Description 

Function  DELAY  uses  a fifth-order  time-dependent  polynomial  to  model  the 
temperature  distribution  in  a pipe  or  duct: 

5 

T(x,t)  = 5 a.(t)x^ 
i=0 


where 


a^  - time-dependent  coefficients  of  the  polynomial 
X = normalized  position  in  the  pipe  or  duct  (0  ^ x < 1) 
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At  time  t'=t+At,  the  temperature  distribution  is  shifted  by  a distance 
Ax=At/r  , where  r„  is  the  transport  time  through  the  conduit: 

X X 

T(x,t')  = T(x-Ax,t)  , Ax  < X i 1 

For  positions  between  x=0  and  x=Ax,  linear  interpolation  is  used  to  evaluate 
t emperatures : 

T(x,t')  = T(0,t')+[T(0,t)-T(0.t')]^  , 0 < X < Ax 

At  each  time  step,  temperatures  are  evaluated  at  five  points  and  a new  set  of 
coefficients  for  the  polynomial  are  found  by  Gaussian  elimination.  It  has 
been  found  empirically  that  for  a stable  solution,  the  evaluation  points 
should  be  equally  spaced  along  the  duct  or  pipe  with  the  following  exceptions: 

1.  One  of  the  points  should  be  located  at  x=Ax. 

2.  If  Ax  > 1/5,  there  should  be  one  or  more  points  in  the  region  0<x<Ax. 

3.  If  Ax  > 1/2,  a linear  temperature  distribution  is  used. 

4.  If  Ax  > 1,  a uniform  temperature  is  used. 

This  method  of  representing  transport  delays  is  necessary  because  the 
simulation  program  selects  its  own  variable  time  step  between  user-supplied 
minimum  and  maximum  values.  Since  a minimum  time  step  of  a tenth  of  a second 
is  not  unreasonable,  representing  transport  delays  by  simply  storing  a table 
of  values  could  require  an  excessively  long  table. 

3.2  Hysteresis  Effects 

Hysteresis  effects  are  modeled  by  the  function  HYSTER.  In  the  context  of 
HVACSIM,  HYSTER  represents  the  relative  position  of  a valve  or  damper,  which 
may  differ  from  the  relative  position  of  the  actuator  controlling  the  valve  or 
damper  position.  Function  HYSTER  has  the  following  arguments: 

C : a control  signal  between  0 and  1,  representing  the  relative  position  of 
an  actuator 

CR  : a variable  related  to  the  position  of  a valve  or  damper 
CRO  : the  value  of  CR  at  the  previous  simulation  time 
TOLD:  the  previous  simulation  time 
HYS  : the  hysteresis  parameter 

C can  be  thought  of  as  the  relative  position  of  an  actuator  which  moves  a 
valve  or  damper,  and  HYS  represents  mechanical  slack  in  the  coupling  between 
the  actuator  and  the  valve  or  damper.  More  specifically,  HYS  is  the  fraction 
of  the  actuator's  range  over  which  HYSTER  remains  constant  when  the  actuator's 
direction  of  travel  reverses.  Like  C,  the  value  of  HYSTER  is  always  between  0 
and  1.  The  relationship  between  C and  HYSTER  is  illustrated  in  figure  3.2.1. 

The  values  of  C and  HYS  are  supplied  by  the  calling  subroutine.  The 
other  three  arguments  should  be  stored  in  the  saved  workspace  by  the  calling 
subroutine.  Thus  the  function  might  be  used  as  in  the  following  line: 

CH  = HYSTER(C,SAVED(1),SAVED(2),SAVED(3),HYS) 
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Figure  3.2.1.  Hysteresis  model. 


3.3  Bessel  Functions,  Polynomial  Curve  Fitting,  and  Fin  Efficiencies 

The  subroutines  described  in  this  section  are  taken,  with  minor  modifica- 
tions, from  reference  [8].  The  first  two  subroutines,  described  brief ly  in 
sections  3.3.1  and  3.3.2,  calculate  the  mod i f ied  Be s se  1 functions  Ij^(x)  and 
Ej^(x).  Section  3.3.3  describes  a subroutine  used  for  polynomial  curve  fit- 
ting. The  subroutine  described  in  section  3.3.4  uses  the  first  three  subrou- 
tines to  find  heat  exchanger  fin  efficiencies  under  any  operating  conditions, 
provided  that  condensation  does  not  occur. 

3.3.1  SUBROUTINE  BESI(X,N,BI,  lER) 


This  subroutine  computes  the  modified  Bessel  function  I^(x).  The  argu- 
ments of  the  subroutine  are  defined  below. 


X:  the  argument  of  the  Bessel  function 

N:  the  order  of  the  Bessel  function  (integer  N 

BI:  the  resultant  value  of  the  Bessel  function 

lER:  the  resultant  error  code: 


lER  = 0 
lER  = 1 
lER  = 2 
lER  = 3 
lER  = 4 


no  error 
N < 0 
X < 0 

BI  < 10“38 
X > 90. 


> 0) 


(BI  set  to  1.) 
(BI  set  to  1.) 
(BI  set  to  0.) 
(BI  set  to  10^8) 
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3.3.2  SUBROUTINE  BESK  (X,  N,BK,  lER) 


This  subroutine  computes  the  modified  Bessel  function  Kj^(x).  The  argu- 
ments of  the  subroutine  are  defined  below. 

X:  the  argument  of  the  Bessel  function 

N:  the  order  of  the  Bessel  function  (integer  N >.  0) 

BK:  the  resultant  value  of  the  Bessel  function 

lER:  the  resultant  error  code: 


3.3.3  SUBROUTINE  POLFIT(X,Y,N,TOL,LAST,COF) 

This  subroutine  fits  a polynomial  to  the  ordered  pairs  of  data  points 
(x,y)  and  returns  the  coefficients  of  the  polynomial.  The  arguments  of  the 
subroutine  are  defined  below. 


TOL:  an  error  tolerance 

LAST:  the  maximum  order  of  the  fitted  polynomial  (1  ^ LAST  ^ 9) 

COF:  an  array  of  (LAST+1)  coefficients  to  the  polynomial 

The  subroutine  begins  by  fitting  a first  order  polynomial  to  the  set  of 
points.  If  the  standard  deviation  of  the  correlation  is  greater  than  the 
error  tolerance  TOL  and  the  order  of  the  polynomial  is  less  than  LAST,  the 
order  is  increased  by  one  and  a new  set  of  coefficients  is  determined.  This 
procedure  is  repeated  until  either  the  error  tolerance  is  met  or  the  maximum 
order  specified  by  LAST  is  reached. 

3.3.4  SUBROUTINE  SUFED(ROH,COF) 

The  purpose  of  this  subroutine  is  to  obtain  coefficients  for  a polynomial 
representing  heat  exchanger  fin  efficiencies  (in  the  absence  of  condensation) 
as  a function  of  the  fin  surface  heat  transfer  coefficient  and  the  fin  geome- 
try. The  first  argument  of  the  subroutine,  ROH,  is  the  ratio  of  tube  outside 
diameter  to  fin  diameter: 

ROH  = Dj„/Dj. 

For  rectangular  plate  fins,  is  the  diameter  of  a circular  fin  having 
the  same  area  per  tube  as  the  rectangular  fin.  The  subroutine  returns  a set 
of  five  coefficients  in  the  array  COF.  These  coefficients  can  be  used  to 
determine  heat  exchanger  dry  fin  efficiencies: 


lER  = 0 
lER  = 1 
lER  = 2 
lER  = 3 
lER  = 4 


N < 0 
X i 0. 

X > 85. 

BI  > 10^^ 


no  error 


(BI  set  to  0.) 
(BI  set  to  0.) 


X:  an  array  of  up  to  60  x-coordinates 

Y:  an  array  of  up  to  60  y-coordinate s 

N:  the  number  of  data  points  (X,Y)  (N  ^ 60) 


^dry  = + ^2F  + CgF^  + C4F^  + CjF^ 
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where 


F = 0.5  (Df  - (2h/k6)°-^ 

h = dry  surface  heat  transfer  coefficient 
k = fin  thermal  conductivity 


6 = fin  thickness 


The  subroutine  determines  the  fin  efficiency 
sixty  operating  points,  using  the  subroutines 
tine  POLFIT  to  determine  the  coefficients. 


by  the  method 
BESI  and  DESK, 


of  Gardner  [9]  at 
and  calls  subrou- 


3.4  Refrigerant  Properti^ 


A number  of  functions  are  available  for  determining  refrigerant  proper- 

tles/C’e"  HoL  aerived  Uo.  T 

cients  for  twelve  different  refrigerants  [lOJ.  me  lasi  aigu 
function  is  an  ind«,  IR,  coding  for  the  refrigerant.  He  meaning 

f ol  lows : 


IR  Refrigerant  IB 


1 

2 

3 

4 

5 

6 


R-11 

R-12 

R-13 

R-14 

R-21 

R-22 


7 

8 

9 

10 

11 

12 


R-23 

R-113 

R-11 4 

R-500 

R-502 

R-318 


He  functions  require  SI  units  on  input,  and  return  values  in  SI  • 
Internally,  inputs  are  converted  to  English  units,  and  results  are  converted 
back  to  SI  units.  The  required  units  are  as  follows: 


Variable 


Symbol  Units 


Pressure 
Specific  Volume 
Temperature 
Entropy 
Enthalpy 
Liquid  Density 


P 

V 

T 

S 

H 

P 


kPa 


m^/kg 


°C 

kJ/(kg»C) 

kJ/kg 

kg/m^ 


Several  of  the  functions  are  iterative  and  require  a relative  error 
tolerance,  RTOL.  With  one  exception,  the  tolerance  is  applied  to  a Fahrenh^t 
temperature.  An  absolute  tolerance  of  0.001  is  built  into  the  function.s.  The 
single  exception  is  the  function  which  returns  specific  volume.  In  this  case, 
the  relative  error  tolerance  RTOL  is  applied  to  the  density  (i.e.,  the  inverse 
of  the  specific  volume),  and  no  absolute  tolerance  is  used.  All  of  the  iter-„ 
ative  functions  are  derived  from  the  basic  equations  in  reference  [10],  using 
the  Newton-Raphe  son  method. 


The  list  of  available  functions  and  their  meanings  is  as  follows: 
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Func  tion 
Name 

( Arguments) 

Value 

Returned 

PSAT 

(T,  IR) 

^sat 

TSAT 

(P,  RTOL,  IR) 

^sat 

TVSAT 

(V,  RTOL,  IR) 

T 

sa  t 

PGAS 

(V,  T,  IR) 

P 

VGAS 

(P,  T,  RTOL,  IR) 

V 

HGAS 

(P,  V,  T,  IR) 

E 

SGAS 

(V,  T,  IR) 

S 

BPS 

(P,  S,  RTOL,  IR) 

H 

TPH 

(P,  H,  RTOL,  IR) 

T 

TVH 

(V,  H,  RTOL,  IR 

T 

DULAT 

(P,  V,  T,  IR) 

^Hiat 

RHOLIQ 

(T,IR) 

P 

CV 

(V,T,IR) 

Cv 

In  addition  there  is  one  subroutine  which  takes  V,  T,  and  IR,  and  determines 
C , C , the  ratio  C„/C„,  and  the  velocity  of  sound  in  the  refrigerant: 

p'  V p V 

SUBROUTINE  CPCV  (V,  T,  IR,  CP,  CV,  GAMMA,  SONIC) 

3 . 5 Steam  and  Liquid  Water  Properties 

The  functions  described  in  this  section  compute  water  and  steam  proper- 
ties. Some  care  is  required  in  their  use,  since  many  of  the  functions  are 
valid  over  a narrow  range  of  conditions  (e.g.  only  at  saturation).  Descrip- 
tions of  the  functions  are  divided  into  four  sections:  superheated  steam, 
saturated  steam,  saturated  liquid  water,  and  water  at  atmospheric  pressure. 
Three  of  the  functions  are  listed  in  both  of  the  last  two  sections.  These 
functions  are  based  on  data  for  saturated  water  at  temperatures  greater  than 
100  °C,  and  water  at  atmospheric  pressure  for  temperatures  from  0 to  100.  The 
latent  heat  of  vaporization  has  been  included  in  the  saturated  steam  category. 

In  all  of  the  water  and  steam  property  functions,  temperatures  are  in  Celsius, 
pressures  are  in  kilopascals,  specific  volumes  are  in  cubic  meters  per  kilo- 
gram, and  entropy  is  in  kilojoules  per  kilogram  per  degree  Celsius. 


3.5.1  Superheated  Steam 


Function  Name 

Value 

and  Arguments 

Returned 

References 

VS 

(P, 

T) 

V 

(m^/kg) 

eq.  from  [11] 

HS 

(P, 

T) 

H 

(kj/kg) 

eq.  from  [11] 

SS 

(P, 

T) 

S 

[kJ/(kg  K)] 

eq.  from  [11] 

CPS 

(T) 

Cp 

[kJ/(kg  K)] 

eq.  from  [12] 

TPSS 

(P, 

S) 

r 

(®C) 

uses  SS  and  CPS 

VISSPH 

(T) 

p 

[kg/(m  s)] 

data  from  [13] 

STEAMK 

(T) 

K 

[W/(m  K)] 

data  from  [13] 

CVS 

(V, 

T) 

Cv 

[kj/(kg  K)] 

adapted  from  [14] 
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3.5.2  Saturated  Steam 


Function  Name 

Value 

and  Arguments 

Returned 

References 

TSATS 

(P) 

T (C) 

eq.  from  [11] 

PS  ATS 

(T) 

P (kPa) 

eq.  from  [11] 

VSATS 

(P,  T) 

V (m^/kg) 

eq.  from  [11] 

HSATS 

(T) 

H (kJ/kg) 

eq.  from  [11] 

SSATS 

(T) 

S [kJ/(kg  K)] 

eq.  from  [11] 

VISSV 

(P) 

p [kg/ (m  s) ] 

data  from  [13] 

HFG 

(T) 

AHfg  (kJ/kg) 

eq.  from  [11] 

3.5.3 

Saturated  Liquid  Water 

Function  Name 

Value 

and  Arguments 

Returned 

References 

VSATW 

(T) 

V (m^/kg) 

eq.  from  [11] 

ESATW 

(T) 

H (kJ/kg) 

eq.  from  [11] 

SSATW 

(T) 

S [kj/(kg  K)] 

eq.  from  [11] 

WMU 

(T) 

p [kg/(m  s)l 

data  from  [15] 

for  T > 

100  C. 

WCP 

(T) 

C [kJ/(kg  K)] 

data  from  [15] 

for  T > 

100  C. 

WK 

(T) 

K [kW/(m  K)] 

data  from  [15] 

for  T > 

100  C. 

3.5.4 

Liquid  Water 

at  1 Atmosphere 

WMU 

(T) 

p [kg/(m  s)] 

eq. 

from 

[16] 

WCP 

(T) 

C [kJ/(kg  K)] 

eq. 

from 

[16] 

WK 

(T) 

K [kW/(m  K)] 

data 

from 

[16] 

WRHO 

(T) 

p (kg/m^) 

eq. 

from 

[16] 

3.6 

Air  Properties 

Five  functions  and  one  subroutine  are  available  for  computing  properties 
of  air  [11],  which  can  be  treated  as  an  ideal  gas.  The  subroutine  argument 
list  is  as  follows: 

SUBROUTINE  CPCVA(TCA,  CPA,  CVA,  GAMMA,  SONIC) 

where  TCA  is  the  Celsius  air  temperature  at  which  the  remaining  arguments  are 
computed,  CPA  and  CVA  are  specific  heats  at  constant  pressure  and  constant 
volume,  respectively,  in  kJ/(kg  K),  GAMMA  is  the  ratio  CPA/CVA,  and  SONIC  is 
the  speed  of  sound  in  meters  per  second. 

The  five  air  property  functions  are  as  follows: 

AKA  (T)  : Thermal  conductivity  [KW/(m  K)  ] as  a function  of  temperature  (C). 

VISCA  (T) : Dynamic  viscosity  [kg/(m  s)]  as  a function  of  temper-ature  (C). 

HA  (T)  : Enthalpy  (kj/kg)  as  a function  of  temperature  (C). 
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PHIA  (T)  The  entropy  function  (defined  below)  [kJ/ (kg  K)]  as  a function  of 
temperature  (C). 

TPHIA  (T):  Temperature  (C)  as  a function  of  the  entropy  function  [kJ/(kg  K)]. 

The  entropy  function,  t,  is  defined  by 
T 

0 = j dT 


where  Tq  is  a reference  temperature  at  which  entropy  is  zero.  Information  on 
the  use  of  the  entropy  function  can  be  found  in  reference  [11]  or  in  standard 
thermodynamics  textbooks. 
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4. 


COMPONENT  MODELS 


4.0  INTRODUCTION 

This  chapter  is  devoted  to  descriptions  of  component  models  developed  for 
use  with  MODSIM.  Wherever  possible,  experimental  or  theoretical  verification 
of  the  models  has  been  included  with  the  descriptions.  In  Table  4.0.1,  the 
available  component  models  are  listed  in  functional  categories.  Several  of 
the  components  (TYPES  12,  13,  16,  and  17)  have  been  included  in  two  categor- 
ies. 

The  remainder  of  this  section  relates  to  the  calculation  of  pressures  and 
flow  rates,  which  is  common  to  many  component  models. 

Pressure  and  Flow  Calculations 


Component  models  generally  relate  pressure  losses  to  flow  rates  using 
equations  of  the  following  form: 

AP  = Kw^ 

where  AP  = pressure  drop  across  the  device 

K = flow  resistance  or  pressure  loss  coefficient 
w = mass  flow  rate 

On  the  basis  of  how  pressure  and  flow  calculations  are  handled,  the  component 
models  can  be  divided  into  four  categories: 

1.  No  fluid  flows  through  the  component. 

2.  Fluid  flows  through  the  component  but  no  pressure  drop  is  calculated. 
Any  flow  resistance  may  be  modeled  by  adding  it  to  the  flow  resis- 
tance of  an  adjacent  component. 

3.  Outlet  pressures  and  inlet  flow  rates  are  used  to  calculate  inlet 

pressures.  If  necessary,  outlet  flow  rates  are  calculated  by  mass 
conservation. 

4.  Inlet  and  outlet  pressures  are  used  to  calculate  flow  rates. 

The  connections  among  pressures  and  flow  rates  in  a simulation  define  a set  of 
equations  which  must  be  solved  simultaneously.  In  setting  up  a simulation, 
the  user  must  be  careful  to  ensure  that  the  set  of  equations  is  neither  under- 
determined nor  overdetermined,  but  has  a unique  solution.  In  general,  one 
component  from  the  fourth  category  is  needed  in  each  flow  circuit,  to  deter- 
mine a flow  rate.  The  inlet  pressure  of  this  component  will  generally  be  a 
boundary  condition  (i.e.  not  an  output  of  any  other  component),  providing  a 
fixed  pressure  basis  for  the  determination  of  all  other  pressures  in  the 
circuit.  For  further  discussion  of  pressure  and  flow  rate  calculations,  refer 
to  the  example  simulation  chapter  of  reference  [1].  In  Table  4.0.2,  the 
component  models  are  organized  into  the  four  categories  listed  above. 
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Table  4.0.1  Component  models  arranged  in  functional  categories 


Fans  and  Pumps 
TYPE  1:  Fan  or  Pxunp 

Conduits 

TYPE  2:  Conduit  (Duct  or  Pipe) 

TYPE  3:  Inlet  Conduit  (Duct  or  Pipe) 

Flow  Splits  and  Merges 
TYPE  4:  Flow  Merge 
TYPE  6:  Flow  Split 

TYPE  13:  Three-way  Valve  with  Actuator 
TYPE  17:  Mixing  Dampers  and  tferge 
TYPE  18:  Plenum 
TYPE  21:  'Grounded'  Flow  Split 

Valves.  Dampers,  and  Flow  Restrictors 

TYPE  5:  Damper  or  Valve 

TYPE  9:  Linear  Valve  with  Actuator 

TYPE  13:  Three-way  Valve  with  Actuator 

TYPE  17:  Mixing  Dampers  and  Merge 

TYPE  23:  Steam  Nozzle 

TYPE  24:  Ideal  Gas  Nozzle 

Sensors  and  Controllers 

TYPE  7:  Temperature  Sensor 

TYPE  8:  Proportional-Integral  Controller 

TYPE  16:  'Sticky'  Proportional  Controller 

TYPE  20:  High  or  Low  Limit  Controller 

TYPE  26:  Control  Signal  Inverter 

Heat  Exchangers 

TYPE  10:  Hot  Water  to  Air  Heating  Coil  (Single) 
TYPE  11:  Hot  Water  to  Air  Heating  Coil  (Detailed) 
TYPE  12:  Cooling  or  Dehumidifying  Coil 
TYPE  25:  Steam  to  Air  Heat  Exchanger 

Humidifiers  and  Dehumidifiers 
TYPE  12:  Cooling  or  Dehumidifying  Coil 
TYPE  14:  Evaporative  Humidifier 
TYPE  22:  Steam  Spray  Humidifier 

Conditioned  Zones 
TYPE  15:  Room 

Superblock  Coupling  Components 

TYPE  16:  'Sticky'  Proportional  Controller 

TYPE  19:  Flow  Balance  Control 
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Table  4.0.2  Component  models  arranged  by  pressure 


/flow  calculation  method 


Category  1:  No  Fluid  Flow 

TYPE  7:  Temperature  Sensor 
TYPE  8:  Proportional-Integral  Controller 
TYPE  16:  'Sticky'  Proportional  Controller 
TYPE  20:  High  or  Low  Limit  Controller 
TYPE  26:  Control  Signal  Inverter 

Category  2:  No  Pressure  Drop 

TYPE  12:  Cooling  or  Dehumidifying  Coil 
TYPE  14:  Evaporative  Humidifier 
TYPE  15:  Room 
TYPE  18:  Plenum 

TYPE  22:  Steam  Spray  Humidifier 

TYPE  25:  Steam  to  Air  Heat  Exchanger  (Air  Side) 

Category  3 : Inlet  Pressure  Calculated  from  Flow  Rate  aud  Outlet  Pressure 

TYPE  1:  Fan  or  Pump 

TYPE  2:  Conduit  (Duct  or  Pipe) 

TYPE  4:  Flow  Merge 

TYPE  5:  Damper  or  Valve 

TYPE  6:  Flow  Split 

TYPE  9:  Linear  Valve  with  Actuator 

TYPE  10:  Hot  Water  to  Air  Heating  Coil  (Simple) 

TYPE  11:  Hot  Water  to  Air  Heating  Coil  (Detailed) 

TYPE  13:  Three-way  Valve  with  Actuator 
TYPE  19:  Flow  Balance  Control 

Category  Flow  Rate  Calculated  from  Inlet  and  Outlgl  Pressures 

TYPE  3:  Inlet  Conduit  (Duct  or  Pipe) 

TYPE  17:  Mixing  Dampers  and  Merge 
TYPE  21:  'Grounded'  Flow  Split 
TYPE  23:  Steam  Nozzle 
TYPE  24:  Ideal  Gas  Nozzle 

TYPE  25:  Steam  to  Air  Heat  Exchanger  (Steam  Side) 


Flow  Resistance  Coefficients 

As  stated  above,  component  models  generally  relate  pressure  losses  to 
flow  rates  using  equations  of  the  following  form: 

AP  = Kw^ 

When  AP  is  in  kPa  and  w is  in  kg/s,  K has  dimensions  of  g'^m'^.  If  AP  is  in 
Pascals,  K is  in  kg’’^m~^. 

A dimensionless  loss  coefficient,  K',  commonly  reported  in  the  litera- 
ture, is  defined  as  follows: 
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head  loss 


2AP 


K'  = 


velocity  head 


pv- 


where  p is  the  density  of  the  fluid  and  v is  its  velocity, 
dimensionless,  AP  must  be  in  pascals  when  p is  in  kg/m^  and  v is 
relationship  between  the  dimensionless  parameter  K'  and  the 
parameter  K is 


For  K'  to  be 
in  m/s.  The 
dimens ioned 


K = 


0.001  K' 
2pA^ 


where  A is  the  cross-sectional  flow  area  and  the  factor  of  0.001  converts  K 
for  use  with  pressures  in  kilopascals. 


K may  also  be  expressed  as  a function  of  the  dimensionless  Darcy  friction 
factor,  f: 


8fL 

K = 

pn^d^ 

where  d is  the  hydraulic  diameter  of  the  pipe  or  duct  and  L is  its  length. 


Valves  and  Dampers 

For  valves  and  dampers,  K is  a function  of  a control  signal  defining  the 
position  of  the  valve  or  damper,  and  of  parameters  defining  the  characteris- 
tics of  the  valve  or  damper.  Two  parameters  common  to  the  valve  and  damper 
models  are  X,  a fractional  leakage  rate  at  unit  AP  when  the  valve  or  damper  is 
closed,  and  a flow  resistance  when  the  valve  or  damper  is  open.  The  flow- 
handling capacity  of  a control  valve  is  commonly  represented  in  the  literature 
by  a flow  coefficient  or  valve  constant,  Cy»  defined  as  follows: 

= (q)(AP)"®-5 

where  q is  the  volumetric  flow  rate  in  gallons  per  minute  and  AP  is  the 
pressure  drop  in  pounds  per  square  inch  when  the  valve  is  fully  open.  When  AP 
is  in  kPa  and  w is  in  kg/s, 

= 1.732  X 10^  (pC^)"2 

A valve  is  essentially  a variable  fluid  resistance.  The  manner  in  which 
the  flow  rate  varies  with  valve  stem  position  is  known  as  the  valve  character- 
istic, or  inherent  characteristic,  since  it  is  always  measured  in  practice 
with  a constant  pressure  drop  across  the  valve.  For  a linear  valve,  the  flow 
rate  is  directly  proportional  to  the  valve  stem  position  when  the  pressure 
drop  across  the  valve  is  constant. 

It  is  important  to  distinguish  between  the  inherent  characteristic  of  a 
valve  and  the  installed  characteristic.  The  latter  depends  on  both  the  in- 
herent characteristic  and  the  authority  of  the  valve,  a: 

, = 

'^y  * *:o  * 
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where 


AP^  = pressure  drop  across  the  open  valve, 

APg  = pressure  drop  across  the  remainder  of  the  circuit, 

IT  = flow  resistance  of  the  open  valve,  and 

Kg  = sum  of  flow  resistances  around  the  remainder  of  the  circuit. 

For  three-way  valves,  APg  and  Kg  refer  only  to  that  portion  of  the  circuit 
where  the  flow  rate  is  being  controlled  by  the  inlet  port  in  question,  and  t e 
two  inlet  ports  may  have  different  authorities. 


The  installed  characteristic  of  a valve  is  strictly  identical  to  the 
inherent  characteristic  only  when  the  authority  is  one,  that  is,  when  the 
valve  represents  the  only  pressure  drop  in  the  system.  A valve  with  low 
authority  has  little  effect  on  the  flow  rate  over  most  of  its  operating  range. 


The  leakage  flow  rate  through  a closed  valve  also  depends  on  the  valve  s 
authority,  as  well  as  on  the  leakage  parameter,  X.  For  the  valve  and  damper 
models  described  in  this  chapter,  the  normalized  or  fractional  leakage  rate  is 
given  by 


^ = )i[a  + !i^(l  - a)]-®-5 
'^max 


Dampers  also  have  inherent  characteristics,  installed  characteristics, 
and  authority.  The  effects  of  valve  and  damper  authority  are  accurately 
represented  by  the  valve  and  damper  models  described  in  this  chapter. 
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4.1  TYPE  1:  FAN  OR  PUMP 


General  Description 

The  pump  or  fan  model  calculates  a pressure  rise  and  an  efficiency  as 
functions  of  a mass  flow  rate,  using  dimensionless  performance  curves.  Hie 
performance  curves  are  represented  by  polynomials,  with  empirically  determined 
coefficients  which  must  be  supplied  by  the  user.  The  efficiency  is  used  to 
compute  a temperature  rise  across  the  fan  or  pump,  as  well  as  power  consump- 
tion. The  last  parameter,  MODE,  determines  whether  air  or  water  properties 
(density  and  specific  heat)  are  to  be  used. 

Nomenclature 


Cf  - dimensionless  flow  coefficient 

Cjj  - dimensionless  pressure  head  coefficient 

Cp  - specific  heat  of  fluid  (kJ/kg-C) 

D*^  - diameter  of  blades  or  impeller  (m) 

E - power  consumption  (kW) 

N - rotational  speed  (rps) 

- inlet  pressure  (kPa) 

Pq  - outlet  pressure  (kPa) 

T^  - inlet  temperature  (C) 

Tq  - outlet  temperature  (C) 
w - mass  flow  rate  (kg/s) 
r\  - efficiency 
p - fluid  density  (kg/m*) 

Mathematical  Description 

The  dimensionless  flow  coefficient  and  the  dimensionless  pressure  head 
coefficient  are  defined  as  follows: 

C.  = 

pNd3 

r _ lOOOAP 
^h 

pN^D^ 

where  AP  is  the  pressure  rise,  P^  - in  kilopascals.  The  factor  of  1000 

converts  the  pressure  difference  from  pascals  to  kilopascals  for  dimensional 
consistency.  For  this  reason,  the  units  given  in  the  nomenclature  section 
must  be  used. 

The  model  determines  the  head  coefficient,  Cj^,  and  the  efficiency,  t),  as 
functions  of  the  flow  coefficient,  C£: 

^ = Sq  + + ^2^f^  ^ ^3^f^  ^ *4^f^ 

+ e^Cf  + e2Cf^  + + e4Cf^ 

jlhe  'a'  and  'e'  coefficients  must  be  supplied  as  parameters  by  the  user.  The 
inlet  pressure  is  then  obtained  from  the  definition  of  Cj^: 
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p.  = 
1 

= Po  - 

O.OOlCiipN^D^ 

The  efficiency 
calculated  on 

is  used  to  compute  power  consumption, 
the  assumption  that  all  inefficiencies 

E = 

wAP 

np 

To  = 

= Ti  + 

(Po  - Pi)  (1  - 1) 

TYPE  1 Component  Conf ieuration 

Inputs 

Description 

1 

w 

- fluid  mass  flow  rate 

2 

Po 

- outlet  pressure 

3 

N 

- rotational  speed 

4 

Ti 

- inlet  temperature 

Outputs 

Description 

1 

Pi 

- inlet  pressure 

2 

To 

- outlet  temperature 

3 

E 

- power  consumption 

Parameters 

Description 

1-5 

- coefficients  for  head  vs.  flow  curve 

6-10 

- coefficients  for  efficiency  curve 

11 

D 

- diameter 

12 

Mode 

- air  = 1 , water  = 2 
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4.2  TYPE  2:  CWJDUIT  (DUCT  OR  PIPE) 


General  Description 

The  conduit  model  is  designed  to  account  for  three  effects:  thermal 
losses  to  ambient  conditions,  transport  delays,  and  dynamics  due  to  thermal 
capacitance.  The  model  has  four  modes,  determined  by  the  last  parameter.  If 
the  absolute  value  of  MODE  is  1,  air  properties  are  used  and  the  model  repre- 
sents a duct.  If  the  absolute  value  of  MODE  is  2,  water  properties  are  used 
and  the  model  represents  a pipe.  The  sign  of  MODE  determines  the  method  used 
to  calculate  thermal  capacitance  effects.  The  merits  of  the  two  methods  are 
discussed  below. 

Nomenclature 


A - surface  area  of  conduit 

Cgj  - thermal  capacitance  (mass  * specific  heat)  of  conduit 
Cp  - specific  heat  of  fluid 
h - heat  transfer  coefficient 
K - flow  resistance  coefficient 
P - pressure 

T - temperature 

t - time 

U - overall  heat  transfer  coefficient 

V - volume  of  conduit 
w - mass  flow  rate  of  fluid 
p - density  of  fluid 
T - time  constant 

- capacitance  term  of  time  constant 

- transport  time 

Sub  scri  2ts 

amb  - ambient 
i - inlet  or  inside 
o - outlet  or  outside 
ss  - steady  state 

Mathematical  Description 

The  pressure  at  the  inlet  of  the  conduit  is  calculated  from  the  outlet 
pressure  and  the  mass  flow  rate,  using  a flow  resistance  coefficient,  K,  which 
is  assumed  constant: 

= Pq  + sign(w)Kw^  (4.2.1) 

Ihe  steady  state  outlet  temperature  is  given  by 

^ss  " T^b  - Tam|,)exp(-y)  (4.2.2) 

where  y = PA 

WLp 
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If  MODE  is  positive,  the  temperature  dynamics  are  represented  by  the  following 
equation: 

flo  ^ IssJlIq  (4.2.3) 

dt  T 


where 


r 


1 _Sn 

hi  + ho  J wCp 


A derivation  of  this  equation  can  be  found  in  reference  [171.  This  equation 
is  solved  analytically  within  the  subroutine,  bypassing  the  MODSIM  non-linear 
differential  equation  solver: 


+ (t;  - Tj^)exp(-  % (4.2.4) 

This  result  is  then  sent  to  the  DELAY  subroutine,  which  returns  the  current 
outlet  temperature: 

Tq  = DELAYCtJ.Tj)  (4.2.5) 

where  is  the  transport  time  through  the  conduit: 


T 


X 


If  MODE  is  negative,  the  temperature  dynamics  are  represented  by  the  following 
equation: 

flo  ^ -exp(-g) { _ DELAYETgg  + Tgexp(-y-o)£Ii] ) (4.2.6) 

dt  T-  dt 

where 


2ho 


and  Y defined  above.  This  equation  is  derived  by  inverse  transformation  of 
an  approximate  transfer  function  given  in  reference  [18],  and  is  solved  out- 
side the  subroutine  by  the  MODSIM  non-linear  differential  equation  solver. 


Model  Verification 


The  approximate  transfer  function  from  which  equation  4.2.6  is  derived 
has  been  compared  with  a more  exact  transfer  function  [18].  The  error  of  the 
approximation  is  a function  of  a and  t^.(o,  where  u)  is  an  oscillation  frequency, 
but  the  maximum  error  is  a function  of  a alone.  When  a is  less  than  or  equal 
to  one,  the  error  never  exceeds  10%  in  magnitude  and  8**  in  phase  angle.  (As 
discussed  below,  some  additional  error  will  be  introduced  by  the  transport 
delay  subroutine.)  For  further  details  see  reference  [18].  Since  a is 
directly  proportional  to  the  length  of  the  conduit,  accuracy  can  be  improved 
by  dividing  a long  conduit  into  two  or  more  shorter  conduits. 
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For  a unit  step  increase  in  the  inlet  temperature,  equation  4.2.6  yields 
an  outlet  temperature  step  increase  of  magnitude  expC—h^A/wC^)  at  the  end  of 
the  delay  time,  followed  by  an  exponential  rise  with  a time  constant  of 
Tpexp(-a).  The  time  constant  of  the  simpler  model,  t,  is  essentially  always 
smaller  than  this,  which  helps  to  compensate  for  the  lack  of  an  output  step 
change  in  the  simpler  model.  A comparison  between  the  two  models  is  provided 
by  figure  4.2.1.  The  parameter  values  in  this  comparison  were  calculated  for 
a 0.508  meter  (20  inch)  square  duct,  65.8  meters  (216  feet)  in  length,  carry- 
ing air  at  1.7  cubic  meters  per  second  (3600  cubic  feet  per  minute).  The 
transport  delay  in  this  example  is  ten  seconds. 


100  200  300 


400  500  600 

TIME  (s) 


700  800  900  1000 


Figure  4.2.1  Conduit  model  response  to  a step  change  in  the  fluid  inlet  tem- 
perature at  time  zero. 
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Limitations  of  the  Model 


Although  a step  change  in  the  inlet  temperature  has  been  applied  here  for 
° f 4-tia  mode  with  changes  in  the  inlet 

purposes  of  toV^  tl»e  fs  uot  generally 

^is^ls  aue  to^t.e^  intrlnsu 

if  rtranaiUd^BodV.  siep  change  In  the  inlet  temperature  ‘ 

in  the  term  sent  to  DELAY,  and  a tamp  input  produces  « st«P  f,' 

The  detailed  mode  is  generally  unreliable 

accuracy  diminishes  when  the  flou  rate  is  such  that  the  transport  ‘t”'  «««* 
the  duration  of  a ramp  input  or  the  period  Jf/“‘hen  tte 

characteristics  of  DELAY  also  limit  the  accuracy  of  the 

transport  time  exceeds  the  time  interval  between  two  step  changes  or  the 
period  of  a sinusoidal  input. 


TYPE  2 Component  Configuration 
Inputs  Description 

w - fluid  mass  flow  rate 


Pq  - outlet  pressure 


- inlet  temperature 


amb 


1 
2 

3 

4 

5 

Outputs 
1 

2 ^i  " inlet  pressure 

Parameters  Description 


- ambient  temperature 

- outlet  temperature 
Description 

- outlet  temperature 


1 

2 

3 

4 

5 

6 


hjA  - inside  surface  heat  transfer  coefficient  ♦ area 

h A - outside  surface  heat  transfer  coefficient  ♦ area 
o 

C - thermal  capacitance  of  conduit  material 
m 

V ~ volume  of  conduit 

K - flow  resistance  coefficient 

H - height  of  outlet  above  inlet 


7 


MODE  - 


-2  =>  water,  detailed  dynamics 
-1  =>  air,  detailed  dynamics 

1 =>  air,  simple  dynamics 

2 =>  water,  simple  dynamics 


(one 

D.E.) 

(one 

D.E.) 

(no 

D.E.) 

(no 

D.E.) 
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4.3  TYPE  3:  INLET  CONDUIT  (DUCT  OR  PIPE) 


General  Description 

In  the  TYPE  2 Conduit  model,  the  inlet  pressure  is  calculated  from  the 
outlet  pressure  and  the  mass  flow  rate.  In  the  TYPE  3 Inlet  Conduit  model, 
the  mass  flow  rate  is  calculated  from  the  inlet  and  outlet  pressures.  In  all 
other  respects  the  two  models  are  identical. 

Nomenclature 


K - flow  resistance  coefficient 
AP  - inlet  pressure  minus  outlet  pressure 
w - mass  flow  rate 
Mathematical  Description 

Equation  4.2.1  in  the  TYPE  2 Conduit  model  is  replaced  by  the  following 
equation : 

w = sign(AP)  (AP/K)®*^ 

Temperatures  are  calculated  as  described  for  the  TYPE  2 Conduit  model. 

TYPE  3.  Component  Configuration 


Inputs 

Description 

1 

Pi 

- inlet  pressure 

2 

Po 

- outlet  pressure 

3 

Ti 

- inlet  temperature 

4 

T 

amb 

- ambient  temperatur 

5 

To 

- outlet  temperature 

Outputs 

Description 

1 

To 

- outlet  temperature 

2 

w 

- mass  flow  rate 
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Parameters 


Description 


1 

2 

3 

4 

5 


lijA  - inside  surface  heat  transfer  coefficient  ♦ area 
hjjA  - outside  surface  heat  transfer  coefficient  ♦ area 
Cjg  - thermal  capacitance  of  conduit  material 
V - volxime  of  conduit 

K - flow  resistance  coefficient 


6 H - height  of  outlet  above  inlet 


7 

MODE  - 

-2  => 

water. 

detailed  dynamics 

(one 

D.E.) 

-1  => 

air. 

detailed  dynamics 

(one 

D.E.) 

1 => 

air. 

simple  dynamics 

(no 

D.E.) 

2 => 

water. 

simple  dynamics 

(no 

D.E.) 

4.4  TYPE  4:  FLOW  MERGE 


Genera  1 Description 

This  component  models  the  merging  of  two  flow  streams.  It  calculates  the 
mass  flow  rate  and  temperature  of  the  single  outlet  flow  stream,  and  the 
pressures  at  the  two  inlets.  The  model  assumes  that  the  flow  resistance 
parameter  is  the  same  for  all  three  branches. 

Nomenclature 


K - flow  resistance  coefficient 
P - pressure 

T - temperature 

w - mass  flow  rate 

Subscripts 

1 - first  inlet 

2 - second  inlet 

3 - outlet 

Mathematical  Description 

The  outlet  mass  flow  rate  is  simply  the  sum  of  the  inlet  flow  rates: 

W3  = Wf  + W2 


The  outlet  temperature  is  a weighted  average  of  the  inlet  temperatures: 

T3  = wiTi  + W2T2 

Inlet  pressures  are  calculated  from  the  outlet  pressure,  the  flow  rates,  and 
the  flow  resistance  coefficient: 

^1  ~ P3  y[sign(w]^) wj^^  + sign(w2)v3^] 

P2  “ ^3  sign(w3)w3^] 

TYPE  4 Component  Configuration 
Description 

- flow  rate  at  first  inlet 

- flow  rate  at  second  inlet 

- pressure  at  outlet 

- temperature  at  first  inlet 

- temperature  at  second  inlet 


Inputs 

1 

2 

3 

4 

5 


^1 

^2 

P3 

Tl 
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Outputs 

Description 

1 W3 

- flow  rate  at  outlet 

2 Pi 

- pressure  at  first  inlet 

3 Pj 

- pressure  at  second  inlet 

4 T3 

- temperature  at  outlet 

Parameters 

Description 

1 K 

- flow  resistance  coefficient 
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4.5  TYPE  5:  DAMPER  OR  VALVE 


General  Description 

A damper  or  valve  is  essentially  a variable  fluid  resistance.  Tbe  manner 
in  which  the  flow  rate  varies  with  damper  or  valve  position  is  known  as  the 
damper  or  valve  characteristic,  or  inherent  characteristic.  The  inherent 
characteristic  is  always  measured  in  practice  with  a constant  pressure  drop 
across  the  damper  or  valve.  For  a linear  valve,  for  example,  the  flow  rate  is 
directly  proportional  to  the  valve  stem  position  when  the  pressure  drop  across 
the  valve  is  constant. 

This  model  represents  dampers  or  valves  having  inherent  characteristics 
which  are  linear,  exponential,  or  intermediate  between  linear  and  exponential. 
The  characteristic  is  determined  by  the  value  of  the  third  parameter,  W^. 
When  W£  is  zero,  the  model  represents  an  exponential  valve  or  damper.  When  Wf 
is  one,  a linear  valve  or  damper  is  modeled.  Intermediate  values  of  Wj  pro- 
duce intermediate  characteristics. 

The  model  computes  the  inlet  pressure  as  a function  of  mass  flow  rate  and 
relative  damper  or  valve  position.  A non-zero  leakage  parameter  is  required 
to  prevent  infinite  flow  resistance  when  the  damper  or  valve  is  fully  closed. 
Parameter  values  for  modeling  opposed  leaf  dampers  and  parallel  leaf  dampers 
are  given  in  the  model  validation  section  below. 

It  is  important  to  distinguish  between  the  inherent  characteristic  of  a 
damper  or  valve  and  the  installed  characteristic.  The  latter  depends  on  both 
the  inherent  characteristic  and  the  authority  of  the  damper  or  valve,  as 
discussed  in  Section  4.0.  The  installed  characteristic  of  an  inherently 
linear  valve  is  strictly  linear  only  when  the  authority  is  one,  that  is,  when 
the  valve  represents  the  only  pressure  drop  in  the  system.  A valve  with  low 
authority  has  little  effect  on  the  flow  rate  over  most  of  its  operating  range. 
When  used  in  a system  with  other  components,  the  TYPE  5 component  accurately 
models  the  effects  of  authority. 

Nomenclature 

C - relative  position 

K - flow  resistance  coefficient  [0.001/ (kg  m) ] 

M - mode:  M = 0 =>  damper  closed  when  C = 0. 

M = 1 =>  damper  closed  when  C = 1. 

P - pressure  (kPa) 

w - mass  flow  rate  (kg/s) 

W£  - weighting  factor  for  linear  term  of  flow  resistance  coefficient 

k - leakage  parameter:  fractional  leakage  when  AP  = 1 kPa. 

Subscripts 

0 - full  open 

1 - inlet 

2 - out let 
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Mathematical  Description 


If  M 0 then  C = 1.  - C 

(2C  - 2) 

J Vo + <1  - w 

[(1  - X)C  + 

Pj  = P2  + sign(vr)Kw^ 

K is  the  flow  resistance  as  a function  of  C,  the  relative  position  (0  < C 
^ 1).  The  inherent  characteristic  of  the  damper  or  valve  model  is  determined 
by  the  parameter  W£.  When  W£  is  zero,  the  model  represents  an  exponential 
valve  or  damper.  When  W£  is  one,  a linear  valve  or  damper  is  modeled.  Inter- 
mediate values  of  W£  produce  intermediate  characteristics. 

Model  Validation 


Reference  [19]  gives  curves  of  dimensionless  flow  resistance  versus 
damper  position  for  parallel  leaf  dampers  and  for  opposed  leaf  dampers,  which 
are  the  two  types  of  dampers  most  commonly  used  in  HVAC  control  applications. 
The  above  equation  for  E was  converted  to  dimensionless  form  as  described  in 
section  4.0,  and  a nonlinear  regression  program  was  used  to  find  values  for 
the  parameters  E^,  X,  and  W£  which  minimized  the  difference  between  the  above 
equation  for  E and  the  published  curves.  The  following  values  were  obtained: 


2pA^E, 


W, 


Opposed 

.52274 

.01091 

.61785 


Parallel 

.54205 

.012214 

.8689 


Reference  [20]  contains  curves  of  flow  versus  position  at  several  authorities 
for  parallel  leaf  and  opposed  leaf  dampers.  As  an  independent  check  on  the 
model,  a simulation  was  run  with  the  inlet  conduit  model  (TYPE  3)  joined  to 
the  TYPE  5 component  with  the  parameters  given  above,  and  the  authority  was 
varied  by  varying  the  conduit  resistance.  Figures  4.5.1  and  4.5.2  compare  the 
simulation  results  with  published  curves  of  the  same  damper  authorities  for 
opposed  and  parallel  leaf  dampers  respectively  [20].  For  parallel  leaf  damp- 
ers, agreement  between  the  model  and  the  published  curves  is  excellent.  For 
opposed  leaf  dampers  the  results  are  less  dramatic,  but  better  agreement  might 
be  obtained  by  decreasing  the  leakage  parameter  used  in  the  simulations. 
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Figure  4.5.1  Flow  versus  position  for  an  opposed  leaf  damper,  as  a function 
of  damper  authority,  A. 


T-W 

Figure  4.5.2  Flow  versus  position  f^r  a parallel  leaf  damper,  as  a function 
of  damper  authority,  A. 
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TYPE  5.  Component  Configuration 


Inputs 

Description 

1 w 

- mass  flow  rate  [ kg/s  ] 

2 

- outlet  pressure  [ KPa  ] 

3 C 

- relative  position  of  damper  or  valve  [ - ] 

Outputs 

Description 

1 Pi 

- inlet  pressure  [ KPa  ] 

Parameters 

Description 

1 ^0 

- full  open  resistance  [ ,001/(kg-m)  ] 

2 X 

- leakage  parameter  [ - ] 

3 Wj 

- weighting  factor  for  linear  term  of  flow  resistance  [ 

4 M 

- mode:  M = 0 =>  damper  closed  when  C = 0. 

M = 1 =>  damper  closed  when  C = 1. 
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4.6  TYPE  6:  FLOW  SPLIT 


General  Description 

This  component  models  the  division  of  a flow  stream  into  two  flow 
streams.  It  calculates  the  mass  flow  rates  at  the  two  outlets,  and  the  pres- 
sure at  the  single  inlet.  The  model  assumes  that  the  flow  resistance  parame- 
ter is  the  same  for  all  three  branches.  Care  has  been  taken  to  ensure  that 
correct  results  are  obtained  even  when  the  flow  direction  of  one  or  more  flow 
streams  is  reversed. 

Nomenclature 


K - flow  resistance  coefficient 
P - pressure 

w - mass  flow  rate 

Subscripts 

0 - center  of  component 

1 - inlet 

2 - first  outlet 

3 - second  outlet 

Mathematical  Description 

The  model  is  derived  from  the  following  set  of  equations: 

2 

P^  - Pq  = 0*5  K sign(wj)  wj 

Pq  - P2  = 0.5  K sign(w2>  ^2^ 

Pq  - P3  = 0.5  K signCw^)  Wj^ 

Wi  = W2  + W3 

where  Pq  is  the  pressure  at  the  center  of  the  component  and  0.5  K is  the  flow 
resistance  between  the  center  and  an  inlet  or  outlet.  The  set  of  equations  is 
solved  for  W2  and  W3,  and  the  solution  is  expressed  in  the  following  form: 

w,  = + [sign(wi)]  ^ if  S > 0. 

2 Kwi 

«2  = + [sign(AP)]  if  s < 0 

2 K 4 


where 

and 


AP  = P3  - 

S = Wf2  - 


2 abs(AP) 


E 


4  5 


W3  = Wf  - W2 


When  S is  positive,  all  three  flow  rates  have  the  same  sign,  and  Wj^  represents 
the  single  inlet  (positive)  or  outlet  (negative)  flow  rate  of  a split  or 
merge.  When  S is  negative,  aiid  Wg  do  not  have  the  same  sign,  and  wj^ 
represents  one  of  two  inlets  (wj  positive)  or  outlets  (wj^  negative). 

The  pressure  at  the  nominal  inlet,  P^,  is  then  calculated: 

Pj  = P3  + y[sign(wj) Wj^  + sign(w3)w3^] 

TYPE  6 Component  Configuration 
Description 

- flow  rate  at  inlet 

- pressure  at  first  outlet 

- pressure  at  second  outlet 
Description 

- flow  rate  at  first  outlet 

- flow  rate  at  second  outlet 

- pressure  at  inlet 
Description 

- flow  resistance  coefficient 


Inputs 
1 
2 
3 

Outputs 
1 
2 
3 

Parameters 


"^1 

P2 

P3 

^^2 

^^3 


K 
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4.7  TYPE  7:  TEMPERATURE  SENSOR 


General  Description 

Temperature  sensors  are  modeled  by  a simple  first-order  differential 
equation  with  a single  time  constant.  The  input  to  the  component  is  modified 
using  a gain  and  an  offset  supplied  as  parameters,  allowing  the  output  to  be 
treated  as  a voltage  level,  a control  signal,  or  a temperature  in  either 
Fahrenheit  or  Celsius.  If  a Celsius  temperature  output  is  desired,  an  offset 
of  zero  and  a gain  of  one  are  used.  If  a control  variable  between  zero  and 
one  is  required  for  use  as  a controller  input,  the  offset  is  the  minimum 
allowable  temperature  and  the  gain  is  the  maximum  allowable  temperature  minus 
the  minimum  temperature. 

Nomenc 1 ature 

C^  - temperature  input  modified  by  gain  and  offset 
Cq  - output  signal  (e.g.  sensed  temperature) 

T^  - temperature  input 

T - gain 

T®  - offset 

At  - time  step 

T - time  constant 


Mathematical  Description 


dt 


C-  - C 
1 o 


If  the  time  constant,  T,  is  greater  than  or  equal  to  one  second,  the  differen- 
tial equation  is  solved  outside  the  subroutine  by  the  MODSIM  differential 

If  r is  less  than  one  second. 


equation  solver, 
solved  within  the  subroutine: 


the  differential  equation  is 


IF  (r/At)  < 0.05  OR 


C = C. 
o 1 


ELSE 


ICi  - c„- 


10  THEN 


Co  = Ci  - (Ci  - )exp(-At/t) 
where  C^  is  the  value  of  C^  one  time  step  ago. 
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TYPE  7 Component  Configuration 


Inputs 

Description 

1 Ti 

- temperature  input 

- temperature  sensor  output 

Outputs 

Description 

- temperature  sensor  output 

Parameters 

Description 

1 X 

- temperature  sensor  time  constant 

2 

- offset 

- gain 
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4.8  TYPE  8:  PROPORTIONAL- INTEGRAL  CONTROLLER 


General  Description 

This  subroutine  models  an  analog  proportional-integral  controller.  Pro- 
portional control  alone  can  be  modeled  by  setting  the  integral  gain  to  zero. 
A time  constant  is  used  to  model  the  response  time  of  the  controller  itself. 
The  feedback  signal,  generally  representing  a temperature,  must  be  a control 
variable,  which  can  be  obtained  by  passing  a temperature  through  the  TYPE  7 
temperature  sensor  model. 

Nomenclature 


- controller  output  (0  ^ C ^ 1) 

- feedback  signal  (controlled  variable) 

- set  point 

- steady  state  controller  output 

- error  signal 

- integral  portion  of  controller  output 

- integral  gain 

- proportional  gain 

- time  constant  associated  with  controller  response 


Mathematical  Description 
The  error  signal  is  defined  by 


^ ^set  ^f 


The  integral  portion  of  the  output  signal  is  given  by 

41 

dt 


4^  = K^E 


The  output  control  signal  is  then  determined: 

dC  C,  - C 
_ ss 

dt  T 


where 


Ccc  = K„E  + I 
ss  p 


This  differential  equation  is  solved  by  MODSIM  unless  the  time  constant,  t,  is 
less  than  one  second,  in  which  case  it  is  solved  within  the  subroutine: 


IF  (x/At)  < 0.05 
C = C, 


OR 


- c I < 10 


-10 


THEN 


ss 


ELSE 


C = Css  - ^^ss  " ^ )exp(-At/x) 
where  C is  the  value  of  C one  time  step  ago. 


TYPE  8 Component  Configuration 


Inputs 

Description 

1 

Cf 

- controlled  variable 

2 

Cset 

- set  point  for  controlled  variable 

3 

I 

- integral  portion  of  control 

signal 

4 

C 

- output  control  signal  (from 

second 

Outputs 

Description 

1 

I 

- integral  portion  of  control 

signal 

2 

C 

- output  control  signal  (0  ^ C 

< 1) 

Parameters 

Description 

1 

•"p 

- proportional  gain 

2 

Ki 

- integral  gain 

3 

V 

- controller  time  constant 

(from  first  output) 
output) 


(0  < I < 1) 
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4.9  TYPE  9:  LINEAR  VALVE  WIIH  PNEUMATIC  ACTUATOR 


General  De scription 

A valve  is  essentially  a variable  fluid  resistance.  Tbe  manner  in  which 
the  flow  rate  varies  with  valve  stem  position  is  known  as  the  valve  character- 
istic, or  inherent  characteristic,  since  it  is  always  measured  in  practice 
with  a constant  pressure  drop  across  the  valve.  For  a linear  valve,  the  flow 
rate  is  directly  proportional  to  the  valve  stem  position  when  the  pressure 
drop  across  the  valve  is  constant. 

It  is  important  to  distinguish  between  the  inherent  characteristic  of  a 
valve  and  the  installed  characteristic.  The  latter  depends  on  both  the  inher- 
ent characteristic  and  the  authority  of  the  valve,  as  discussed  in  Section 
4.0.  Tbe  installed  characteristic  of  an  inherently  linear  valve  is  strictly 
linear  only  when  the  authority  is  one,  that  is,  when  the  valve  represents  the 
only  pressure  drop  in  the  system.  A valve  with  low  authority  has  little 
effect  on  the  flow  rate  over  most  of  its  operating  range.  When  used  in  a 
system  with  other  component  models,  the  TYPE  9 valve  accurately  models  the 
effects  of  valve  authority.  Figure  4.9.1  illustrates  the  effect  of  valve 
authority  on  the  installed  characteristic  of  a linear  valve.  This  figure  was 
obtained  using  the  TYPE  9 linear  valve  model  in  conjunction  with  a TYPE  3 pipe 
model . 


Figure  4.9.1  Linear  valve  characteristics  as  a function  of  authority,  A. 
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Hie  valve  model  includes  a pneumatic  actuator  model.  Actuator  hysteresis 
is  represented  using  the  utility  subroutine  HYSTER<  described  in  section  3.2. 
The  dynamic  response  of  the  actuator  is  modeled  by  a first  order  differential 
equation.  The  actuator  can  be  effectively  removed  by  setting  the  time  con- 
stant and  the  hysteresis  parameter  to  zero. 

When  the  flow  rate  through  the  valve  goes  to  zero,  the  pressure  drop 
across  the  valve  becomes  indeterminate.  To  prevent  this  problem,  a non-zero 
leakage  parameter  is  required  by  the  model. 

Nomenclature 

C - requested  relative  valve  position  (0  ^ C ^ 1) 

- relative  actuator  position  (0  j<  <.  1) 

Cjj  - relative  valve  position  (0  ^ Cj^  ^ 1) 

hys  - fraction  of  actuator's  range  over  which  remains  constant 
when  actuator's  direction  of  travel  reverses 
K - flow  resistance  parameter  when  valve  is  open  (Cj^  = 1) 
w - mass  flow  rate 

X - leakage  parameter:  fractional  leakage  when  AP  = 1. 

X - actuator  time  constant 

Mathematical  Description 

The  relationship  between  the  input  control  signal,  C,  and  the  actuator  posi- 
tion, Cg,  is  defined  by: 


dt  T 

This  differential  equation  is  solved  by  MODSIM  unless  the  time  constant,  t,  is 
less  than  one  second,  in  which  case  the  following  solution  is  used: 

IF  (r/At)  < 0.05  OR  - C"|  < 10“^®  IHEN 

Ca  = C 


ELSE 


= C - (C  - C^“)exp(-At/T) 
where  is  the  value  of  one  time  step  ago. 

The  actuator  position  differs  from  the  valve  position,  Cj^,  due  to  hysteresis 
effects,  which  are  determined  by  the  utility  function  HYSTER  described  in 
sect  ion  3.2 . 

Cj^  = HYSTER(C^,  hys) 

The  inlet  pressure  is  given  by 

= Pq  + sign(w)Kw^[(l-k)Cj^  + 
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TYPE  2.  Component  Conf ignration 


Inputs 

Description 

1 Po 

- pressure  at  outlet 

2 w 

- mass  flow  rate 

3 C 

- input  control  signal  (0  <.  C ^ 1) 

- actuator  relative  position 

Outputs 

Description 

- actuator  relative  position 

2 P, 

- pressure  at  inlet 

Parameters 

Description 

1 K 

- flow  resistance  when  valve  is  open 

2 T 

- actuator  time  constant 

3 X 

- leakage  parameter 

4 hy  s 

- hysteresis  parameter 
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4.10  TYPE  10;  HOT  WATER  TO  AIR  HEATING  COIL  (SIMPLE) 

General  Description 

The  TYPE  10  subroutine  is  the  simplest  of  three  water-to-air  heat  exchan- 
ger models  available  at  present  in  the  TYPES  library.  A constant  overall  heat 
transfer  coefficient  (DA)  is  assumed.  Steady  state  air  and  water  outlet 
temperatures  are  calculated  using  an  approximate  effectiveness  relationship 
for  cross-flow  heat  exchangers.  Dynamic  outlet  temperatures  are  calculated 
using  first-order  differential  equations,  in  which  the  time  constant  is  a 
function  of  both  the  coil  flush  time  and  the  thermal  capacitance  of  the  coil. 
However,  the  outlet  temperature  is  not  delayed  by  the  flush  time  of  the  coil. 

Nomenclature 


C — capacitance  rate  (specific  heat  times  mass  flow  rate)  of  air 
Cj^in  - minimum  of  and 

C - capacitance  rate  (specific  heat  times  mass  flow  rate)  of  water 
K - flow  resistance  parameter  on  air  side  of  coil 

- flow  resistance  parameter  on  water  side  of  coil 
NTD  - number  of  transfer  units 

R — ratio  of  minimum  to  maximum  capacitance  rate 

UA  - overall  heat  transfer  coefficient 

w„  - mass  flow  rate  of  air 
& 

w^  - mass  flow  rate  of  water 

e - effectiveness 

T - coil  time  constant 

Tg  - capacitive  term  of  coil  time  constant 

- coil  flush  time 

Subscripts 

a - air 

i - inlet 

o - outlet 

ss  - steady  state 

w - water 


Mathematical  Description 

The  heating  coil  model  uses  an  approximate  equation  for  effectiveness  as 
a function  of  NTD  (number  of  transfer  units)  for  a cross-flow  heat  exchanger 
with  both  fluids  unmixed  [21]  to  determine  steady  state  air  and  water  outlet 
temperatures : 

e = 1 - exp  { [exp(-Rn{NTD})-ll /Rn} 

where 

n = (NTD)"®*2^ 
and  NTD  is  defined  by 
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NTU  = 


UA 


min 


The  steady  state  air  and  water  outlet  temperatures  are  found  using  the  defini- 
tion of  e : 


T = T • + (T  • - T -)eC  • /C 
aoss  ai  '^wi  ai^  min'  a 

T = T . - fT  - T .)C  /C 

woss  wi  ' aoss  ai'  a'  w 

A coil  time  constant  is  also  calculated  by  the  model,  using  the  following 
equa  tion : 

T = ( 

where  is  a capacitive  term  supplied  as  a parameter  and  is  the  coil  flush 
time,  which  is  determined  from  the  coil  volume: 

. = fiVol 


w, 


w 


The  equation  for  the  time  constant  was  chosen  to  give  the  correct  qualitative 
behavior  of  T as  a function  of  water  flow  rate  (see  reference  [22],  figure  7). 
Air  and  water  dynamic  outlet  temperatures  are  both  calculated  using  this  time 
constant : 


dT. 


wo  _ 


dt 


dT 


ao  _ 


T - T 
woss  wo 


T - T 
aoss  ao 


dt 


If  both  and  the  coil  volume  are  zero,  the  differential  equations  are 
eliminated  from  the  model  and  the  steady  state  temperatures  are  used. 

Finally,  inlet  pressures  are  calculated  from  the  outlet  pressures  and  flow 
rates : 


Po4 

= Po,^ 

+ K„w 

ai 

ao 

a i 

P . 

Wl 

= P 
wo 

+ V, 

w 
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TYPE  10  Component  Configuration 


Inputs 

Description 

1 Pao 

- outlet  air  pressure 

- air  mass  flow  rate 

3 T,i 

- inlet  air  temperature 

4 T^o 

- outlet  air  temperature  (from  first  output) 

5 Pwo 

- outlet  water  pressure 

^ ** 

- water  mass  flow  rate 

7 T . 

Wl 

- inlet  water  temperature 

8 T 

wo 

- outlet  water  temperature  (from  second  output) 

Outputs 

Description 

1 

ao 

- outlet  air  temperature 

2 T^o 

- outlet  water  temperature 

3 Pai 

- inlet  air  pressure 

4 P • 

Wl 

- inlet  water  pressure 

Parameters 

Description 

1 UA 

- overall  heat  transfer  coefficient  (times  area) 

2 

- flow  resistance  parameter,  air  side 

3 ^ 

- flow  resistance  parameter,  water  side 

4 Vol 

- water  side  volume  of  coil 

3 

- capacitive  term  of  coil  time  constant 
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4.11  TYPE  11:  HOT  WATER  TO  AIR  HEATING  COIL  (DETAILED) 


General  Description 

Like  the  TYPE  10  component  model,  this  subroutine  represents  a water-to- 
air  cross-flow  heating  coil.  Unlike  TYPE  10,  a constant  overall  heat  transfer 
coefficient  is  not  assumed:  separate  inside  and  outside  heat  transfer  coeffi- 
cients are  calculated  as  functions  of  the  water  and  air  flow  rates.  In 
addition,  a more  complicated  method  is  used  to  represent  temperature  dynamics. 
The  differential  equations  used  to  calculate  temperature  dynamics  in  this 
model  are  derived  from  transfer  functions  for  the  responses  of  the  air  and 
water  outlet  temperatures  to  perturbations  in  the  water  inlet  temperature. 


Nomenclature 


NTU  - 
R 


s - 
UA  - 


w 

w 

e 


a 

w 


T 


Subscripts 


a 

i 

o 

ss 

w 


capacitance  rate  (specific  heat  times  mass  flow  rate)  of  air 
thermal  capacitance  (specific  heat  times  mass)  of  coil  material 
minimum  of  C„  and  C„ 

capacitance  rate  (specific  heat  times  mass  flow  rate)  of  water 
flow  resistance  parameter  on  air  side  of  coil 
flow  resistance  parameter  on  water  side  of  coil 
number  of  transfer  units 

ratio  of  minimum  to  maximum  capacitance  rate 
Laplace  operator 

overall  heat  transfer  coefficient 

mass  flow  rate  of  air 

mass  flow  rate  of  water 

ef  fectiveness 

fin  efficiency 

coil  time  constant 

capacitive  term  of  coil  time  constant 
coil  flush  time 


air 

inlet 

outlet 

steady  state 
water 


Mathematical  Description 


Separate  inside  and  outside  heat  transfer  coefficients  are  calculated  as 
functions  of  the  water  and  air  mass  flow  rates,  w^  and  w„: 

W a 


(qhA)^  = aw^^ 
(hA)^  = cw^*^ 


The  parameters  a,  b,  c,  and  d may  be 
calculated  from  existing  correlations 
which  are  generally  given  in  terms  of 
transfer  coefficient  is  then  given  by 


determined  from  experimental  data  or 
(e.g  reference  [13],  pp.  341  and  351), 
the  Reynolds  number.  The  overall  heat 
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UA  = [ (TihA)^'^  + OiA)^"^  ]“^ 

and  the  NTU  (number  of  transfer  units)  value  is,  by  definition, 

UA 

NTU  = 

^in 

If  a constant  UA  is  desired,  b and  d may  be  set  to  zero. 

The  coil  model  uses  an  approximate  equation  for  effectiveness,  e,  as  a 
function  of  NTU,  for  a cross-flow  heat  exchanger  with  both  fluids  unmixed 
[21]: 


e = 1 - exp{  [exp(-Rn{NTU})-l] /Rn} 
where  n = (NTU)“®*^^ 

The  steady  state  air  and  water  outlet  temperatures  are  found  using  the  defini- 
tion of  e : 


T = T . + (T  • - T .)eC  • /C 
aoss  ai  ' wi  ai'  min'  a 


T = T . - (T  - T .)C  /C 
woss  wi  ' aoss  ai'  a'  w 


Coil  time  constants  are  also  calculated  by  the  model,  using  an  assumption 
that  all  dynamics  are  due  to  changes  in  the  water  inlet  temperature.  Gartner 
and  Harrison  [23]  give  the  following  approximate  transfer  function  for  the 
water  outlet  temperature  response  to  perturbations  in  the  water  inlet  tempera- 
ture: 

- ^wo^*^  Cexp(-y)  + exp(-Pj^)T^s]exp(-T^s) 

‘ 

where 


^ exp(  ) 

a 2a 


« = Pa  + ^^2 

2+^4 

Y = Pi  - bh. 

a 

Pi 

Cw 

P2  = 

Cm 

P3  = ^^^w^x 


m 


p ^ (tihA) 


a 
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and  Ti,  C^,  and  are  defined  in  the  nomenclature  section.  and  ^ 

represent  perturbations  from  some  base  temperature.  Using  the  air  inlet 
temperature  as  the  base  temperature,  they  are  replaced  by  (T^q  ~ and 

- Tg^),  respectively.  Treating  ^ as  a constant,  and  replacing  the  steady 
state  solution  with  the  steady  state  temperature  calculated  from  the  coil 
effectiveness,  the  following  equation  can  be  derived  by  inverse  transformation 
of  the  transfer  function: 


dT, 


wo 


dt 


= <Tai  - + DELAY [T„«p(-Pi)flwi  + 

^w  dt 


where  DELAY  is  the  transport  delay  function  described  in  section  3.1.  Under 
some  conditions  the  accuracy  and  reliability  of  the  model  may  be  limited  by 
the  characteristics  of  DELAY.  For  further  discussion  of  the  limitations  see 
section  4 .2. 


A similar  procedure  is  applied  to  the  air  side  of  the  coil.  Bhargava  et 
al.  [24]  give  an  approximate  transfer  function  of  the  following  form: 


«2 


= = 

T;i(s) 
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l+a2T^s+a3T^‘'s 


Reference  [24]  gives  closed-form  expressions  for  the  constants  a2»  and  a^, 
which  are  derived  from  a more  exact  transfer  function  by  the  Fade  polynomial 
method.  This  approximate  transfer  function  can  be  inverse  transformed  into  a 
second  order  differential  equation,  which  can  be  separated  into  a pair  of 
coupled  first  order  equations.  Assuming  a constant  air  inlet  temperature  as 
before,  the  steady-state  air  temperature  can  be  substituted  into  the  equa- 
tions. The  resulting  pair  of  first-order  differential  equations  is 


dt 


(T  - T ) 
^ aoss  ao^ 


dTao  _ 
dt 


^2^ao 


^3^a 


where 


&2  = 


o + Pi 
ay 


y ’ e 


1 


= a 


2 _ (a  + Pji^) 


2 _ 


ay 


+ e,  [ y 


2 2 

ay 


+ X’ 
Y 


+ 2) 


-i  ] 


= 1 + 


Pl?3 


_ exp(-y) 
l-exp(-y) 


and  the  parameters  a,  y,  and  the  P parameters  are  defined  above.  T^  is  a 
dummy  variable,  introduced  in  the  process  of  separating  the  second  order 
equation  into  two  first  order  equations. 
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Model  Validation 


The  TYPE  11  heating  coil  model  was  tested  by  using  experimentally  meas- 
ured temperatures  and  flow  rates  as  inputs  to  the  model,  and  comparing  the 
measured  air  and  water  outlet  temperatures  with  the  calculated  values.  In 
figure  4.11.1,  simulated  and  measured  air  and  water  outlet  temperatures  are 
compared  for  an  experiment  in  which  step  changes  in  the  control  valve  position 
were  used  to  produce  sharp  changes  in  the  water  inlet  temperature.  In  figure 
4.11.2,  the  dynamics  are  due  to  step  changes  in  the  water  flow  rate.  In  all 
cases  the  agreement  between  simulated  and  measured  outlet  temperatures  is 
quite  good.  Similar  results  are  obtained  for  step  changes  in  the  air  inlet 
temperature  and  the  air  flow  rate.  At  least  some  of  the  differences  must  be 
attributed  to  experimental  error,  since  the  rate  of  energy  transfer  calculated 
from  the  measured  water  flow  rate  and  steady  state  temperatures  is  not  quite 
equal  to  the  energy  transfer  rate  calculated  from  the  measured  air  flow  rate 
and  steady  state  temperatures. 


0 2 4 6 8 10  12  14  16 

TIME  (minutes) 

Figure  4.11.1  Heating  coil  response  to  rapid  changes  in  the  water  inlet  tem- 
perature . 
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TEMPERATURE  (“O 


Figure  4.11.2  Heating  coil  response  to  rapid  changes  in  the  water  flow  rate. 


TYPE  11  Component  Configuration 
Inputs  Description 


1 

2 

3 

4 

5 

6 

7 

8 

9 


- outlet  air  pressure 

a U 


- air  mass  flow  rate 


- inlet  air  temperature 

Tao  “ outlet  air  temperature  (from  first  output) 


P^o  “ outlet  water  pressure 


w 


- water  mass  flow  rate 


T^i  - inlet  water  temperature 

Two  ~ outlet  water  temperature  (from  second  output) 

- dummy  temperature  used  to  calculate  air  temperature  dynamics 
(from  third  output) 
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Outputs 

Description 

- outlet  air  temperature 

- outlet  water  temperature 

3 

- dummy  temperature  used  internally  to  calculate  air  tempera* 
ture  dynamics 

4 Pai 

- inlet  air  pressure 

5 P • 

Wl 

- inlet  water  pressure 

Parameters 

Description 

1 a 

- coefficient  for  the  beat  transfer  correlation  (nbA).  = aw„” 

2 b 

- coefficient  for  the  heat  transfer  correlation  (i)hA).  = aw_^ 

3 c 

- coefficient  for  the  heat  transfer  correlation  (bA)^  = 

4 d 

- coefficient  for  the  heat  transfer  correlation  (bA)^  = ®^w^ 

5 

- flow  resistance  parameter,  air  side 

- flow  resistance  parameter,  water  side 

7 Vol 

- water  side  volume  of  coil 

00 

- mass  times  specific  beat  of  coil  material 
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4.12  TTPE  12:  COOLING  OR  DEHTJMIDIFYING  COIL 


General  Description 

This  model  represents  a circular  finned  or  continuous  finned  serpentine 
heat  exchanger  with  four  or  more  rows  in  counterflow  crossflow  configuration. 
The  subroutine  is  based  on  a program  written  by  Elmahdy  and  Mitalas  [8],  with 
minor  modifications  to  convert  it  into  a TYPE  subroutine  and  to  add  dynamics 
to  their  steady  state  model. 

The  subroutine  uses  equations  for  log  mean  temperature  difference  (LMID) 
and  log  mean  enthalpy  difference  (LMHD)  which  are  strictly  correct  only  for 
counterflow  heat  exchangers.  When  four  or  more  rows  are  arranged  in  crossflow 
counterflow  configuration,  treating  the  coil  as  a counterflow  heat  exchanger 
is  a good  approximation.  Use  of  the  model  to  represent  coils  with  fewer  than 
four  rows  is  not  recommended. 

The  model  accounts  for  condensation  on  the  outside  surface  of  the  coil. 
There  are  three  possible  conditions  for  the  coil:  all  wet,  partially  wet,  or 
all  dry.  The  subroutine  determines  which  of  these  conditions  applies,  and 
treats  each  case  separately. 

TYPE  12  requires  a relatively  large  number  of  parameters  describing  the 
physical  characteristics  of  the  coil.  Fin  efficiencies,  heat  transfer  coeffi- 
cients, and  time  constants  are  calculated  within  the  subroutine  from  this 
information.  The  units  given  in  the  nomenclature  section  ^st  be  used. 

Nomenclature 


m 

Di 

G 

h 

H 

J 

Pr 

Q 

Re 

w 

U 

Uc 

V 

w 

T1 

T 

(I) 


- minimum  area  for  air  flow  (m'^) 

- external  surface  area  (m^) 

- slope  of  enthalpy  vs.  temperature  line  at  saturation,  AH/AT 
(kJ  kg"lc"2) 

- correction  factor  (dimensionless) 

- specific  heat  (kJ  kg~^C~^) 

- total  thermal  capacitance  of  coil  material  (kJ/C) 

- tube  inside  diameter  (m) 

- maximum  air  mass  flux,  (kg  m”^s“^) 

- heat  transfer  coefficient  (kJ  kg~^C~^) 

- enthalpy  (kJ  kg“^) 

- Colburn  heat  transfer  J-factor  (dimensionless) 

- Prandtl  number  (dimensionless) 

- heat  transfer  rate  (kW) 

- Reynolds  number  (dimensionless) 

- bulk  water  temperature 

- overall  heat  transfer  coefficient,  temperature  basis  (kW  m 

- overall  heat  transfer  coefficient,  enthalpy  basis  (kg  m 

- velocity  (m/s) 

- mass  flow  rate  (kg/s) 

- fin  efficiency  (dimensionless) 

- time  constant  (s) 

- humidity  ratio  (kg  water  / kg  dry  air) 


C"l) 
2 3-1; 
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Subscripts 

a - air 

d - dry 

i - inside,  inlet 

o - outside,  outlet 

w - water,  wet 


Math^atical  Description 

Tbe  mathematical  model  has  been  described  in  detail  elsewhere  [8.25].  A 
brief  summary  is  provided  below. 

The  average  heat  transfer  coefficient  on  the  water  side,  h^,  is  evaluated 
using  the  turbulent  flow  relation  [26]: 


h^  = 1.429  [1+0.0146  T^] 

with  all  quantities  in  the  units  indicated  in  the  nomenclature  section.  Under 
dry  coil  conditions,  the  average  heat  transfer  coefficient  on  the  air  side, 
hgj,  is  found  from  a correlation  of  the  following  form: 


Jf  “ Re 


'ad 


= J 


G „ Pr 
p,  a 


-2/3 


The  quantities  and  C2  are^  constants  for  a particular  coil,  and  are  calcu- 
lated as  functions  of  the  coil  geometry  [27].  Under  wet  coil  conditions,  the 
dry  coil  heat  transfer  coefficient  is  modified  as  follows: 


h = Cuh  j 
aw  f ad 

where  the  correction  factor  C£  is  a function  of  the  Reynolds  number  [28]. 

Separate  correlations  are  used  to  determine  the  dry  fin  efficiency, 
[9],  and  the  wet  fin  efficiency,  [29].  In  addition  to  the  inside  and 
outside  heat  transfer  coefficients  and  the  fin  efficiency,  the  overall  heat 
transfer  coefficient  includes  terms  representing  the  thermal  resistance  of  the 
tubes,  and  a constant  fouling  factor  of  5 x 10”*  ®C*in^/W  (3  x 10”^ 
®F-ft^/Btu) . 

Under  dry  coil  conditions,  the  steady  state  air  and  water  outlet  tempera- 
tures are  found  by  simultaneous  solution  of  the  following  three  equations: 


® '^w^p,w  ^^o  ” \i^ 

Q = UA.  (LMTD) 
o 


where 
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LMTD 


~ 

When  the  coil  is  wet,  enthalpies  are  used  in  place  of  temperatures: 


Q = % 

Q = NEeJ.  (H,„  - H^i) 

b 


Q = U (LMHD) 

V O 

where 


LMHD 


^^ai  ~ ^o^  ~ ^^ao  ~ ^i^ 

- H.o>'-  l“<=ao  - Hwi) 


H^^  and  are  saturated  air  enthalpies  calculated  at  the  water  temperatures 

Twi  and  T^q»  respectively,  according  to  the  following  relation: 


H = H + b T 
w o w 


The  quantities  and  b are  determined  iteratively. 

Initially  the  entire  air  side  surface  of  the  coil  is  assumed  to  be  wet. 
Based  on  this  assumption,  the  outlet  air  and  water  temperatures  and  the  air 
side  surface  temperatures  at  inlet  and  outlet  are  calculated.  If  the  surface 
temperature  at  the  inlet  is  lower  than  the  inlet  air  dew  point  temperature, 
the  coil  is  in  fact  all  wet.  If  the  surface  temperature  at  the  air  outlet  is 
higher  than  the  inlet  air  dew  point  temperature,  the  coil  surface  is  all  dry. 
If  neither  of  these  conditions  is  met,  the  surface  is  partly  wet  and  partly 
dry.  In  this  case  an  iterative  procedure  is  used  to  find  the  position  in  the 
air  flow  direction  at  which  the  surface  temperature  is  equal  to  the  dew  point 
temperature.  This  position  is  the  boundary  between  the  dry  and  wet  sections 
of  the  coil. 

Dynamics  have  been  added  to  the  steady  state  model  in  a very  simple  and 
somewhat  artificial  manner: 


dT'  T„^  - T!  ^ 
ao  _ ao ao 

dt  T 


dt  T 


ao  _ 
dt 


(I) 


- 


ao 


ao 


where 


X 


Note  that  the  dynamic  variables  T^^,  T^^,  and  are  used  only  in  the  three 
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differential  equations  given  above.  Steady  state  variables  are  used  in  deter- 
mining whether  the  coil  is  wet  or  dry. 

Model  Validation 

The  steady  state  cooling  coil  model  has  been  tested  extensively  by  its 
authors.  Comparisons  between  the  steady  state  model  and  experimental  data  for 
a variety  of  coils  and  operating  conditions  can  be  found  elsewhere  [30,31]. 

The  dynamics  of  the  coil  model  have  been  investigated  experimentally  for 
changes  in  the  water  flow  rate  and  for  changes  in  the  water  inlet  temperature. 
In  both  experiments,  measured  air  and  water  inlet  temperatures  and  a measured 
water  flow  rate  were  used  as  inputs  to  the  model.  The  air  flow  rate  was 
calculated  from  the  measurements  under  steady  state  conditions.  In  both  cases 
the  outside  surface  of  the  coil  remained  dry.  Measured  outlet  temperatures 
are  compared  with  calculated  temperatures  in  figures  4.12.1  and  4.12.2. 
Agreement  between  the  model  and  the  experiments  is  excellent. 


TIME  (seconds) 


Figure  4.12.1  Cooling  coil  response  to  changes  in  the  water  flow  rate. 
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TIME  (seconds) 


Figure  4.12.2  Cooling  coil  response  to  changes  in  the  water  inlet  temperature 


XyPE  12  Component  Configuration 
Inputs  Description 

- water  mass  flow  rate 

- inlet  water  temperature 

- dry  air  mass  flow  rate 

- inlet  air  dry  bulb  temperature 


w 


■^wi 


“^ai 


- inlet  air  humidity  ratio 

“ dynamic  outlet  water  temperature 
T^o  ” dynamic  outlet  air  temperature 

- dynamic  outlet  air  humidity  ratio 
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Outputs 


Description 


1 

2 

3 

4 

5 

6 


T» 

wo 


T» 

^ao 


(i> 


ao 


- dynamic  outlet  water  temperature  [C] 

- dynamic  outlet  air  temperature  [C] 

- dynamic  outlet  air  humidity  ratio  [ - 1 

- total  cooling  load  at  steady  state  [kW] 

- sensible  cooling  load  at  steady  state  [kW] 

- wet  fraction  of  total  surface  area  [ - ] 


Parameters  Description 


1 

2 

3 

4 

5 

6 

7 

8 
9 

10 

11 

12 

13 

14 

15 

16 

17 

18 


c 

A. 


- coil  type;  0 = flat  continuous  fins,  1 = circular  fins 

- primary  (tube  exterior)  surface  area  [m^] 

- secondary  (fin)  surface  area  [m^l 

- internal  surface  area  [m**] 


N. 


N. 


N, 


A^/Af  - ratio  of  minimum  flow  area  to  face  area  [ - ] 
kf  - thermal  conductivity  of  fin  material  [kW  m~^C~^] 

A£  - coil  face  area  [m^] 

- number  of  fins  per  centimeter 

- number  of  tubes  per  row 

- number  of  rows 

- tube  outside  diameter  [m] 

- tube  inside  diameter  [m] 

- fin  thickness  [m] 

- mass  times  specific  heat  of  coil  material  [kJ/C] 

- tube  longitudinal  spacing  (in  air  flow  direction)  [m] 

- fin  diameter  (if  I^.  = 1)  or  fin  length  (if  = 0)  [m] 

- coil  depth  in  air  flow  direction  [m] 

- thermal  conductivity  of  tube  material  [kW  m~^C“^] 


5 

Si 

D,. 
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4.13  TYPE  13:  IHREE-WAY  VALVE  WITH  ACTUATOR 


General  Description 

This  component  represents  a valve  with  two  inlet  ports  and  one  outlet 
port.  The  position  of  the  valve  is  determined  by  an  input  control  signal,  C. 
Port  1 of  the  valve  is  closed  when  C=0  and  open  when  C=l.  The  sixth  parame- 
ter, MODE,  determines  the  inherent  characteristics  of  the  two  inlet  ports.  If 
MODE=0,  both  ports  are  linear.  If  M0DE=1,  port  1 is  exponential  (equal  per- 
centage) and  port  2 is  linear.  If  M0DE=2,  both  ports  are  exponential.  As 
discussed  in  Section  4.0,  the  actual  behavior  of  the  valve  model  in  a system 
simulation  will  depend  on  the  valve's  authority  over  the  system. 


The  valve  model  includes  a pneumatic  actuator  model.  Actuator  hysteresis 
is  represented  using  the  utility  subroutine  HYSTER,  described  in  chapter  3. 
The  dynamic  response  of  the  actuator  is  modeled  by  a first  order  differential 
equation.  The  actuator  can  be  effectively  removed  by  setting  the  time  con- 
stant and  the  hysteresis  parameter  to  zero. 

When  the  flow  rate  through  the  valve  goes  to  zero,  the  pressure  drop 
across  the  valve  becomes  indeterminate.  To  prevent  this  problem,  a non-zero 
leakage  parameter  is  required  by  the  model. 

Nomenclature 


C 

Ca 

^hl 
Cii2 
hy  s 

- 

^2  - 
MODE 


’"1 

’^2 

X 

X 


- requested  relative  valve  position  (0  < C < 1) 

- relative  actuator  position  (0  ^ < 1) 

- relative  opening  of  port  1 (0  ^ Cj^  £ 1) 

- relative  opening  of  port  2 (0  ^ < 1) 

- fraction  of  actuator's  range  over  which  remains  constant  when 

actuator's  direction  of  travel  reverses 

- flow  resistance  parameter  for  inlet  port  1,  when  port  1 is  open 

(Ch  = 1) 

- flow  resistance  parameter  for  inlet  port  2,  when  port  2 is  open 
(Ci  =•  0) 

determines  port  1 / port  2 inherent  characteristics: 

0 =>  linear/linear 

1 =>  exponential/linear 

2 =>  exponential/exponential 
mass  flow  rate  at  inlet  port  1 
mass  flow  rate  at  inlet  port  2 
mass  flow  rate  at  outlet  port  3 

leakage  parameter:  fractional  leakage  when  AP  = 1. 
actuator  time  constant 


Mathematical  Description 

The  relationship  between  the  input  control  signal,  C,  and  the  actuator  posi- 
tion, C«,  is  defined  by: 

B 

f£a  , C - 

dt  T 


This  differential  equation  is  solved  by  MODSIM  unless  the  time  constant,  t,  is 


less  than  one  second,  in  which  case  the  following  solution  is  used: 

IF  (x/At)  < 0.05  OR  ICgg  - C"l  < 10*"^®  IREN 

C = C 
''a 


ELSE 


C-  = C - (C  - C.  )exp(-At/t) 


where  is  the  value  of  one  time  step  ago. 

The  actuator  position  differs  from  the  valve  position  due  to  hysteresis  ef- 
fects, which  are  determined  by  the  utility  function  HYSTER  described  in  sec- 
tion 3.2. 

Cj^  = HYSTER(Cg,  hys) 

Ch2  = 1-  “ CjjL 

The  outlet  flow  rate  is  the  sum  of  the  inlet  flow  rates,  and  the  outlet  tem- 
perature is  a weighted  average  of  the  inlet  temperatures: 

^^3  " ^1  '”2 

"^3  = <^1^1  + ^2^2^  / ^3 

Calculation  of  the  inlet  pressures  depends  upon  the  mode. 

Mode  = 0: 

= P3  + sign(wj)KjWj^[  (l-X)Cjj^  + X]~^ 

Pj  = P3  + sign(w2)K2W2^[(l-X)Cjj2  + 

Mode  = 1: 


-2C 

Pj  = P3  + sign(wj^)KjW2^X  ^ 

P2  = P3  + sign(w2>K2W2^[  (1~X)Cjj2  + X] 
Mode  = 2: 


-2 


-2C 

Pj  = P3  + sign(Wj^)KjWj^X  ^ 

-2C 

P2  = P3  + 8ign(w2>K2W2^X 
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ryPE  13  Component  Configuration 


Inputs 

Description 

1 Wj 

- mass  flow  rate,  inlet  port  1 

2 wj 

- mass  flow  rate,  inlet  port  2 

3 P3 

- pressure  at  outlet 

4 T, 

- temperature,  inlet  port  1 

5 Tj 

- temperature,  inlet  port  2 

6 C 

- input  control  signal  (0  <.  C < 1) 

- actuator  relative  position 

Outputs 

Description 

- actuator  relative  position 

2 W3 

- outlet  mass  flow  rate 

3 Pi 

- pressure  at  inlet  port  1 

4 Pj 

- pressure  at  inlet  port  2 

5 T3 

- temperature  at  outlet 

Parameters 

Description 

1 

- port  1 flow  resistance  when  port  1 is  open 

2 K2 

- port  2 flow  resistance  when  port  2 is  open 

3 X 

- leakage  parameter 

4 T 

- actuator  time  constant 

5 hys 

- hysteresis  parameter 

6 Mode  - determines  characteristics  of  ports  1 and  2 

0 =>  linear/linear 

1 =>  exponential/ linear 

2 =>  exponential/exponential 
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4.14  TYPE  14:  EVAPORATIVE  HUMIDIFIER 


General  Description 

This  component,  which  is  based  on  a model  developed  by  Chi  [32],  repre- 
sents an  air  stream  passing  over  a water  pan  or  through  a porons  pad,  increas- 
ing the  humidity  and  decreasing  the  temperature  of  the  air.  The  model  assumes 
that  the  process  is  adiabatic,  which  means  that  the  water  temperature  is  as- 
sumed to  be  equal  to  the  entering  air  wet  bulb  temperature.  In  addition,  the 
process  is  assumed  to  occur  at  constant  pressure.  Any  pressure  drop  through 
the  humidifier  may  be  modeled  by  increasing  the  flow  resistance  of  an  adjacent 
component. 

The  evaporative  humidifier  model  requires  the  inlet  air  temperature,  hu- 
midity ratio,  and  mass  flow  rate  as  inputs.  It  computes  the  exit  humidity 
ratio,  temperature,  and  mass  flow  rate.  The  wet  bulb  temperature  and  total 
enthalpy  of  the  air  remain  constant  while  the  moisture  content  rises  and  the 
exit  dry  bulb  temperature  decreases.  The  mass  flow  rate  increases  slightly 
due  to  the  addition  of  moisture  to  the  air. 

Nomenclature 


C - specific  heat  at  constant  pressure 
- enthalpy 

hA  - water  to  air  heat  transfer  coefficient  times  area 
T - temperature 

w - mass  flow  rate 

0)  - humidity  ratio 

Subscripts 

a - air 

fg  - latent  heat  of  vaporization 

i - inlet 

o - outlet 

ref  - reference  point  at  which  H = 0 (-20® C) 
s - saturation 

V - water  vapor 

wb  - wet  bulb 


Mathematical  Description 

The  total  enthalpy  at  the  inlet  is  determined  by  the  following  equation  : 

H.l  = - ^r.f)  * ^ - Tr.f> 

1 + “ai 

Next,  the  inlet  wet  bulb  temperature  is  calculated,  using  a correlation  from 
reference  [8]: 

T^b  " 30-9185  - 39.682  C + 20.5841  - 1.758 

where  C = logg(H^^) 
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The  outlet  saturation  humidity  ratio  is  calculated  next,  using  the  subprogram 
PSATS  described  in  Chapter  3 to  find  the  saturation  pressure  at  the  inlet  wet 
bulb  temperature.  The  constant  in  the  following  equation  is  the  water  to  air 
molecular  weights  ratio. 

0.62198  PSATS(T^Ij) 

“ P.tl.  - PSAISd,^) 


where  is  101.325  kPa.  The  outlet  air  humidity  ratio,  temperature,  and 

mass  flow  rate  can  now  be  determined. 


%o  = ‘"ai  <“aos  " *^ai^<l  " ® 


where  R = 


hA 

pa  ai 


The  outlet  air  temperature  and  mass  flow  rate  aregiven  by 

<i)„^  - (I). 


T = T • - 
ao  ai 


ao  ai 
^pa 


AHf 


g 


1+0) 


'ao 


ao 


1 + 0)  • 
ai 


TYPE  14  Component  Configuration 


Inputs 
1 
2 
3 

Outputs 
1 
2 
3 

Parameters 

1 


0) 


ai 


ai 


Description 

- inUt  air  temper.tura  I C ] 

- inlet  air  humidity  ratio  [ - ] 

- inlet  moist  air  mass  flow  rate  [ kg/s  ] 
Description 

- outlet  air  humidity  ratio  [ - ] 

- outlet  air  temperature  [ C ] 

- outlet  moist  air  mass  flow  rate  [ kg/s  ] 
Description 


0) 


ao 


ao 


w 


ao 


hA  ~ water  to  air  heat  transfer  coefficient  times  area 
[ KW/®C  ] 
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4.15  TYPE  15:  ROOM  MODEL 


General  Description 

The  room  model  has  been  described  in  reference  [33],  and  is  loosely  based 
on  reference  [34].  The  model  considers  the  air  mass,  an  interior  mass,  and  a 
wall  mass.  No  modeling  of  the  building  shell  is  attempted,  and  heat  flows  are 
treated  as  inputs  to  the  model. 

The  air  mass  is  artificially  divided  into  two  parts:  a fully  mixed  part 
near  the  ventilation  supply,  and  a piston  flow  part  near  the  ventilation  ex- 
haust. The  fraction  of  tbe  air  mass  in  these  two  parts  is  adjustable.  Each 
part  of  tbe  air  mass  exchanges  heat  with  two  hinds  of  room  mass:  an  interior 
mass  and  a wall  mass.  The  temperature  of  each  room  mass  is  assumed  to  be  uni- 
form. Two  heat  flows  are  also  defined:  a heat  flow  into  the  wall  mass, 
representing  conduction  through  the  walls,  and  a direct  heat  flow  into  the 
fully  mixed  air  mass,  representing  internal  gains. 


Nomenclature 


- specific  heat  of  air 

- specific  heat  of  interior  mass 

- specific  heat  of  wall  mass 

- fraction  of  air  mass  which  is  fully  mixed 

- heat  transfer  coefficient  times  area  of  interior  mass 

- heat  transfer  coefficient  times  area  of  wall  mass 

- interior  mass 

- wall  mass 

- heat  flow  due  to  internal  gains 

- conduction  heat  flow  into  wall  mass 

- ventilation  air  inlet  temperature 

- exhaust  air  temperature 

- interior  mass  temperature 

- temperature  of  fully  mixed  portion  of  air  mass 

- temperature  of  piston  flow  portion  of  air  mass 

- average  room  air  temperature 

- wall  mass  temperature 

- volume  of  room  air  mass 

- mass  flow  rate  of  ventilation  air 

- density  of  air 


Mathematical  Description 


The  fully  mixed  air  mass  temperature  is  given  by  the  following  equation: 
dT 

(fpvc,)  —a  = Qj  + - V + fH,(T,  - i;)  + fHidj  - i^) 

dt 


The  piston  flow  part  of  the  air  mass  exchanges  heat  with  the  interior  mass  and 
the  wall  mass  as  it  moves  through  the  room  according  to  the  following  equa- 
tion: 


= (1  - f)[H^(T^  - Tp(x))  + Hjdj  - Tp(x))] 


wLC^ 


dTp(x) 

dx 


where  L is  the  total  distance  moved  by  tbe  air  mass.  The  solution  to  tbis 
equation  is: 

Tj,(x)  = T,  - (T3  - Vexp(-ax/L) 


where 


a 


(1  - f) 


and  Tp(x)  is  the  temperature  of  the  air  as  a function  of  the  distance  trav- 
eled. The  spatial  average  air  mass  temperature  is  approximately  given  by: 

dt  pV(l  - f) 


where 


T 


s 


- V 


(1  - exp(-a) ) 
a 


The  room  mass  temperatures  are  given  by  the  following  equations: 

M.C,  = 0,  + H,(Tr  - T.) 
dt 


MjCj  fll  = Hj(Tj.  - Tj) 
dt 


Tr  = fT^  + (l-f)Tp 

The  exhaust  air  temperature  is  equal  to  the  temperature  of  the  piston  flow 
portion  of  the  air  mass  at  x=L: 

T.o  = T,  - (T,  - T„)exp(-a) 

A transport  delay  is  applied  to  the  exhaust  air  temperature  using  the  utility 
subroutine  DELAY  with  a delay  time  equal  to  the  flush  time  for  the  piston  flow 
fraction  of  the  room  volume. 


TYPE  IS  Component  Conf ignration 


Inputs 

Description 

1 w 

- mass  flow  rate  of  ventilation  air 

2 T^i 

- ventilation  air  inlet  temperature 

3 \ 

- temperature  of  fully  mixed  portion  of  air  mass 

- wall  mass  temperature 

5 Tj 

- interior  mass  temperature 

« \ 

~ spatial  average  temperature  of  the  piston  flow  portion  of  the 
air  mass 

7 % 

- conduction  heat  flow  into  wall  mass 

00 

© 

M 

- heat  flow  due  to  internal  gains 

Outputs 

Description 

- temperature  of  fully  mixed  portion  of  air  mass 

2 

- wall  mass  temperature 

3 Tj 

- interior  mass  temperature 

4 

- spatial  average  temperature  of  the  piston  flow  portion  of  the 
air  mass 

5 T, 

- average  rocm  air  temperature 

6 ’^ao 

- exhaust  air  temperature 

Parameters 

Description 

1 V 

- volume  of  room  air  mass 

2 ~ thermal  capacitance  of  walls 

3 ~ thermal  capacitance  of  interior  mass 


4 «w 

- heat  transfer  coefficient  times  area  for  the  wall  mass 

5 Hj 

- heat  transfer  coefficient  times  area  for  the  interior  mass 

6 f 

- fraction  of  the  air  mass  which  is  fully  mixed 
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4.16  TYPE  16:  'STICKY'  PROPORTIONAL  CONTROLLER 


General  Description 

This  component  represents  a proportional  controller  with  upper  and  lower 
limits  on  the  rate  of  change  of  the  output  control  signal.  The  limits  are 
determined  by  the  second  parameter,  B.  If  the  error  signal  changes  by  an 
amount  less  than  B,  the  change  is  ignored  and  the  output  control  signal 
remains  constant  until  the  cumulative  changes  in  error  equal  or  exceed  B.  A 
change  in  the  error  signal  with  a magnitude  greater  than  2B  is  treated  as  a 
change  in  error  of  magnitude  B.  Ulus  for  purposes  of  calculating  the  output 
control  signal,  all  changes  in  error  lie  between  B and  2B.  The  measured 
input  signal,  generally  representing  a temperature,  must  be  a control  variable 
between  zero  and  one,  which  can  be  obtained  by  passing  a temperature  through 
the  TYPE  7 temperature  sensor  model. 

NOTE:  this  component  is  to  be  used  only  to  pass  control  signals  between 
superblocks.  The  variable  defined  by  its  output  should  not  be  solved  simulta- 
neously. 

Nomenclature 


B - error  band 

C - output  control  signal  (0  C ^ 1) 

Cthj  “ measured  temperature  input  control  signal  (0  <.  < 1) 

Cjg  - temperature  setpoint  control  signal  (0  < Cj.g  < 1) 

E - effective  error  at  current  simulation  time 

E^  - effective  error  at  previous  simulation  time 

Kq  - offset 

Kp  - proportional  gain 

Mathematical  Description 


AE  - Cjg  Ejj 


If  (IaeI  < B) 

If  (B  < IaeI  < 2b) 
If  (2b  < IaeI) 


E - 

E = Eq  + (B)sign(AE) 


C = Ko  + KpE 

If  (C  < 0)  C = 0 
If  (C  > 1)  C = 1 
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TYPE  16  Component  Conf ianration 


Inputs 

Description 

1 ^Tm 

- measured  temperature  input  control  signal 

2 Cts 

- temperature  setpoint  control  signal 

Outputs 

Description 

1 C 

- controller  output 

Parameters 

Description 

- proportional  gain 

2 B 

- error  band 

3 Ko 

- offset 
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4.17  TYPE  17:  MIXING  DAMPERS  AND  MERGE 


General  Description 

TYPE  17  represents  a pair  of  mixing  dampers,  mechanically  linked  so  that 
one  damper  opens  as  the  other  closes.  The  first  inlet  is  closed  when  the 
input  control  signal,  C,  is  zero,  and  open  when  C is  one.  The  reverse  is  true 
for  the  second  inlet.  Hysteresis  effects  and  an  actuator  time  constant  are 
not  included  in  the  model. 

Flow  resistances  are  calculated  in  the  same  manner  as  in  the  TYPE  5 
damper  or  valve  model.  For  a discussion  of  damper  characteristics  and  recom- 
mended parameter  values,  see  section  4.5. 

The  TYPE  17  component  takes  three  pressures  as  inputs  and  calculates 
three  flow  rates.  In  this  respect  it  is  similar  to  the  TYPE  21  flow  split  and 
the  TYPE  3 inlet  conduit,  both  of  which  also  compute  an  inlet  flow  rate  rather 
than  an  inlet  pressure. 

Nomenclature 


C - control  signal  (0  ^ C ^ 1) 

Kq  ~ flow  resistance  coefficient  when  damper  is  open 
P - pressure 

R - flow  resistance  as  a function  of  damper  position 
T - temperature 

w - mass  flow  rate 

W£  - weighting  factor  for  linear  term  of  flow  resistance 
Subscripts 

1 - first  inlet 

2 - second  inlet 

3 - outlet 


Mathematical  Description 

R^  and  R2  are  the  flow  resistances  of  the  two  inlets: 


R 


1 = 


R2  = 


^f^o 


-2(1-0 


[(1  - X)C  + X]2 

WfKp 

[(1  - X)(l  - C)  + X]2 


+ (1  - 


-2C 


(1  - 


These  equations  model  linear  damper  characteristics  when  W£  is  one, 
exponential  characteristics  when  W£  is  zero,  and  intermediate  characteristics 
for  intermediate  values  of  W£.  These  flow  resistances  are  used  to  find  two  of 
the  three  flow  rates: 
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= sign(Pj  ” P3)  i?l ?3-!-  j 

v?2  = sign(P2  - P3)  [ j 

h 

The  third  flow  rate  is  the  sum  of  the  first  two: 

W3  = wi  + W2 

The  outlet  temperature  is  a weighted  average  of  the  two  inlet  temperatures 


T3 

= »lTl 

+ WjTj 

^3 

TYPE  17 

Component  Conf isuration 

Inputs 

Description 

1 

Pi 

- pressure  at  first  inlet 

2 

P2 

- pressure  at  second  inlet 

3 

P3 

- pressure  at  outlet 

4 

Tl 

- temperature  at  first  inlet 

5 

T2 

- temperature  at  second  inlet 

6 

C 

- input  control  signal 

Outputs 

Description 

1 

"^1 

- flow  rate  at  first  inlet 

2 

»2 

- flow  rate  at  second  inlet 

3 

*3 

- flow  rate  at  outlet 

4 

T3 

- temperature  at  outlet 

Parameters 

Description 

1 Kjj  - full  open  resistance  [ ,001/(kg-m)  ] 

2 X - leakage  parameter  [ - ] 

3 ^f  ~ weighting  factor  for  linear  term  of  flow  resistance  [ - ] 
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4.18  TYPE  18:  PLENUM 


General  Description 

The  output  of  this  component  is  simply  the  sum  of  its  inputs.  The  single 
parameter,  n,  specifies  the  number  of  input  flow  rates  to  be  summed.  Up  to 
ten  flow  rates  may  be  added  together.  Ten  input  flow  rates  are  required  by 
the  model,  but  only  the  first  n are  summed.  The  remaining  inputs  are  ignored. 

Nomenclature 

n - number  of  input  mass  flow  rates  (1  ^ n 10) 

W£  - i^^  input  mass  flow  rate  (kg/s) 

w^^^  - output  mass  flow  rate 

Mathematical  Description 


*out  = *1  + ^2  + . . . + Wj^ 

NOTE:  Inputs  (n  + 1)  through  10  must  be  specified  in  the  simulation,  but  are 

ignored.  Only  the  first  n inputs  are  used. 

TYPE  18  Component  Configuration 

Inputs  Description 

1 ^1  ~ first  mass  flow  rate 

2 ^2  ~ mass  flow  rate 


n ~ mass  flow  rate 


10  Wj^Q  - tenth  mass  flow  rate 


Outputs 

1 


Description 

Wq^^  - outlet  mass  flow  rate 


Parameters 
1 n 


Description 

- number  of  input  mass  flow  rates 
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4.19  TYPE  19:  FLOW  BALANCE  CONTROL 


General  Description 

This  subroutine  does  not  represent  any  physical  component.  It  serves 
instead  to  minimize  mathematical  problems  which  arise  when  a closed  loop  flow 
stream  is  divided  between  two  superblocks.  It  functions  by  providing  one 
superblock  with  an  estimate  of  the  net  flow  resistance  encountered  by  the  flow 
stream  in  another  superblock.  The  purpose  and  use  of  TYPE  19  is  most  easily 
explained  by  way  of  an  example. 

Figure  4.19.1  represents  a very  simple  flow  circuit  consisting  of  two 
ducts,  a fan,  and  a damper.  Figure  4.19.2  shows  one  possible  way  of  modeling 
the  system,  with  half  of  the  circuit  in  one  block  and  half  in  another.  The 
arrows  indicate  the  flow  of  information  among  units  and  blocks.  For  the 
purpose  of  this  example,  only  pressures  and  flow  rates  are  shown.  Pressure  P^ 
is  a boundary  condition,  providing  a baseline  from  which  pressure  changes  can 
be  calculated.  Pressure  P2  and  flow  rate  wj  are  determined  simultaneously 
within  block  1.  Pressure  P^  is  solved  within  block  2.  Pressure  P3  is  deter- 
mined by  simultaneous  solution  of  blocks  1 and  2,  provided  that  the  two  blocks 
are  in  the  same  superblock.  However,  MODSIM  assumes  that  couplings  between 
superblocks  are  weak  enough  so  that  each  superblock  can  be  treated  as  an 
independent  subsystem.  No  simultaneous  equations  are  defined  between  super- 
blocks. Thus  if  blocks  1 and  2 are  in  different  superblocks,  they  will  be 
solved  sequentially  rather  than  simultaneously  and  a self-consistent  solution 
to  the  full  set  of  equations  defining  pressures  and  flow  rates  is  unlikely  to 
be  achieved. 


Fan 


Damper  Flow 
Resistance 


Figure  4.19.1  Schematic  diagram  of  a simple  closed-loop  system. 
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Figure  4.19.2.  Information  flow  diagram  for  a simulation  of  the  system  shown 
in  figure  4.19.1. 


In  Figure  4.19.3,  the  same  system  is  modeled  using  the  Type  19  flow 
balance  component.  Flow  rates  Wj^  and  i>2  are  calculated  essentially  independ- 
ently in  the  two  blocks,  although  both  represent  the  same  flow  stream.  Type 
19  maintains  internally  a variable  flow  resistance  representing  the  net  flow 
resistance  in  block  2.  At  each  time  step  it  examines  the  difference  between 
w^  and  W2,  aiid  modifies  its  internal  resistance  in  such  a way  as  to  reduce  the 
difference  to  zero.  It  then  uses  this  resistance  to  calculate  P3 , which  it 
ultimately  involved  in  determining  both  w.  and  W2«  Note  that  in  both  Figure 
4.19.2  and  Figure  4.19.3,  the  same  • amount  of  information  passes  between  super- 
blocks:  a flow  rate  is  exchanged  for  a pressure.  The  important  difference  it 
that  in  the  latter  case,  the  Type  19  component  actively  takes  responsibility 
for  bringing  the  two  superblocks  to  convergence. 
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P3 


Figure  4.19.3.  Information  flow  diagram  for  a simulation  of  the  same  system 

in  two  superblocks,  using  the  TYPE  19  component. 


Nomenclature 


K - flow  resistance  parameter  at  current  simulation  time 

Kp  - flow  resistance  parameter  at  previous  simulation  time 

- upstream  pressure,  where  flow  stream  exits  this  SB 
Pq  - downstream  pressure,  where  flow  stream  enters  this  SB 
w - mass  flow  rate  determined  in  this  SB 

w'  - mass  flow  rate  determined  in  another  SB 

AP  - upstream  pressure  minus  downstream  pressure 

Aw  - w'  minus  w 

SB  - superblock 
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Mathematical  Description 


By  definition. 


K 


from  which 

dK  ^ -2AP  = z2K 
dw  ^3  w 

The  change  in  K which  will  eliminate  the  difference  between  w'  and  w is  ap- 
proximated by 

AK  = 

w 


or 


AK  = 2K(1  - 1^) 
and  the  new  value  of  K is 
K = Kp  + AK 

The  upstream  pressure  can  then  be  calculated: 

= Pjj  + Kw^ 


TYPE  19  Component  Configuration 
Inputs  Description 

w - mass  flow  rate  determined  in  this  superblock 

p^  - downstream  pressure,  where  flow  stream  enters  this  superblock 
w'  - mass  flow  rate  determined  in  another  superblock 


Outputs 


1 


Pi 


Description 

- upstream  pressure,  where  flow  stream  exits  this  superblock 


Parameters 


1 


K, 


Description 

- flow  resistance  parameter  at  time  = 0. 
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4.20  TYPE  20:  HIGH  LIMIT  OR  LOW  LIMIT  CONTROLLER 


General  Description 

This  component  represents  a proportional  controller  which  exerts  a con- 
trol action  only  when  the  measured  input  signal  is  on  one  side  of  a fixed 
setpoint.  For  a high  limit  controller,  no  control  action  is  exerted  unless 
the  measured  input  signal  exceeds  the  setpoint.  For  a low  limit  controller, 
control  action  is  exerted  only  when  the  controlled  input  falls  below  the  set- 
point.  The  setpoint  is  determined  by  the  absolute  value  of  the  second  parame- 
ter. The  component  represents  a high  limit  controller  when  the  second  parame- 
ter is  positive,  and  a low  limit  controller  when  the  second  parameter  is 
negative.  The  measured  input  signal,  generally  representing  a temperature, 
must  be  a control  variable  between  zero  and  one,  which  can  be  obtained  by 
passing  a temperature  through  the  TYPE  7 temperature  sensor  model.  The 
setpoint  parameter  must  also  be  a normalized  control  variable  with  an  absolute 
value  between  zero  and  one. 

Nomenclature 


C - output  control  signal  (0  ^ C 1 1) 

^Tm  “ measured  temperature  input  control  signal  (0<.Cj.jjj<l) 
Cj.g  - temperature  setpoint  parameter  (0  < < 1) 

Ep  - proportional  gain 

Mathematical  Description 

E = ICjgl  - Cjm 

High  limit  controller  2.  0)* 

It  (Ct„  > Cj.,)  C = C + Kj,E 

Low  Limit  Controller  (Cj.^  < 0): 

If  < \Cj^\)  C = C + KpE 

TYPE  20  Component  Configuration 


Inputs 

2 C 

Outputs 
1 C 

Parameters 


Description 

- measured  temperature  input  control  signal 

- controller  output 
Description 

- controller  output 
Description 

- proportional  gain 

- setpoint  times  i,  where  i = -1  for  low  limit  controller,  or 

i = 1 for  high  limit  controller 
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4.21  TYPE  21:  'GROUNDED'  WATER  OR  AIR  SPLIT 


General  Description 

This  component  models  the  division  of  a flow  stream  into  two  flowstreams. 
It  is  referred  to  as  'grounded'  by  analogy  to  a grounded  electrical  junction: 
the  pressure  at  the  center  of  the  split  is  a constant,  entered  as  a parameter. 
Using  as  inputs  the  pressures  at  the  two  outlets,  the  subroutine  calculates 
the  pressure  and  mass  flow  rate  at  the  single  inlet  and  the  mass  flow  rates  at 
the  two  outlets.  The  model  assumes  that  the  flow  resistance  parameter  is  the 
same  for  all  three  branches.  Care  has  been  taken  to  ensure  that  correct 
results  are  obtained  even  when  the  flow  direction  of  one  or  more  flow  streams 
is  reversed. 

Nomenclature 

K - flow  resistance  coefficient 

P - pressure 

w - mass  flow  rate 


Subscripts 


0 

1 

2 

3 


center  of  split 
inlet 

first  outlet 
second  outlet 


Mathematical  Description 


0.5 


E 


0.5 


W3  = sign(P(j  - P3) 


wi  = W2  + W3 
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TYPE  21  Component  Configuration 


Inputs 

Description 

1 

- pressure  at  first  outlet 

2 P3 

- pressure  at  second  outlet 

Outputs 

Description 

1 Pi 

- pressure  at  inlet 

2 

- flow  rate  at  inlet 

3 ^2 

- flow  rate  at  first  outlet 

4 W3 

- flow  rate  at  second  outlet 

Parameters 

Description 

1 Po 

- pressure  at  center  of  split 

2 K 

- flow  resistance  coefficient 
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4.22  ITPE  22  : STEAM  SPRAY  HUMIDIFIER 


General  Description 

The  steam  spray  humidifier  model  calculates  an  outlet  air  temperature, 
flow  rate,  and  humidity  ratio  for  a constant  pressure  process  in  which  steam 
is  injected  into  an  air  flow  stream  to  increase  the  humidity  of  the  air.  The 
outlet  air  humidity  ratio  is  limited  either  by  the  steam  flow  rate  or  by  the 
saturation  humidity  ratio  (at  atmospheric  pressure)  times  a saturation  effi- 
ciency. The  mass  flow  rate  of  the  air  increases  slightly  due  to  the  addition 
of  moisture. 


Nomenclature 


C 

cP^ 

^pax 

^a  tm 

T • 
ax 

T 

ao 


a X 
*^ao 


a X 
^ao 


aos 


- specific  heat  of  dry  air  [kJ/(kg  K)] 

- specific  heat  of  moist  inlet  air  [kj/(kg  K)] 

- specific  heat  of  steam  or  water  vapor  [kJ/(kg  K)  ] 

- atmospheric  pressure  (101.325  kPa) 

- temperature  of  inlet  air  (C) 

- temperature  of  outlet  air  (C) 

- temperature  of  inlet  steam  (C) 

- mass  flow  rate  of  moist  inlet  air  (kg/s) 

- mass  flow  rate  of  moist  outlet  air  (kg/s) 

- mass  flow  rate  of  inlet  steam  (kg/s) 

- saturation  efficiency  ( - ) 

- inlet  air  humidity  ratio  (kg  water  / kg  dry  air) 

- outlet  air  humidity  ratio  (kg  water  / kg  dry  air) 

- saturation  humidity  ratio  (kg  water  / kg  dry  air) 


Mathemat ical  Description 


pax 


where  C 


The  inlet  air  specific  heat  depends  on  the  inlet  air  humidity  ratio: 
C, 


= ^pa  **^ai^ps 

1 + (i)  . 

ax 


-pg  is  the  specific  heat  of  dry  air  and  Cp^  is  the  specific  heat  of 
steam  or  water  vapor.  The  air  outlet  temperature  is  a weighted  average  of  the 
inlet  air  and  steam  temperatures: 


= W,,C„.,Te  + 


Lao 


s’^ps  s 


ax^pax  ax 


^sSs  """aiSai 


Assuming  that  all  of  the  entering  steam  is  absorbed  by  the  air,  the  outlet 
humidity  ratio  is 


%o  = “ai 


w 


ax 


The  saturation  humidity  ratio  at  the  outlet  temperature  is  calculated  next, 
using  the  subprogram  PSATS  described  in  Chapter  3 to  find  the  saturation 
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pressure  at  the  outlet  temperature.  The  constant  in  the  following  equation  is 
the  ratio  of  water  to  air  molecular  weights. 


0.62198  PSATS(Tao) 

- PSATS(T 
a tm  ao 

where  P^tm  101.325  kPa.  The  outlet  humidity  ratio.  Wgo*  cannot  exceed  the 
saturation  efficiency  times  the  outlet  humidity  ratio  at  saturation,  ilWaoj: 

“ao  = “ao'  i1“aos> 

TYPE  22  Component  Configuration 


Inputs 

Description 

1 

T 

•^s 

- steam  temperature  ( C ] 

2 

T . 
ai 

- inlet  air  temperature  [ C ] 

3 

*^ai 

- inlet  air  humidity  ratio  [ - ] 

4 

^ai 

- inlet  air  mass  flow  rate  [ kg/s  ] 

5 

'^s 

- steam  mass  flow  rate  [ kg/s  ] 

Outputs 

Description 

1 

T 

ao 

- outlet  air  temperature  [ C ] 

2 

'*^ao 

- outlet  air  mass  flow  rate  [ kg/s  ] 

3 

‘^ao 

- outlet  air  humidity  ratio  [ - ] 

Parameters 

1 

- saturation  efficiency 
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4.23  TYPE  23  : STEAM  NOZZLE 


General  Description 


The  steam  nozzle  calculates  a mas 
given  the  downstream  pressure  and  the 
and  pressure.  The  model  assumes  an 
condensation  into  account  and  uses  r< 
steam  properties. 


s flow  rate  and  exit  steam  temperature, 
upstream  steam  stagnation  temperature 
adiabatic,  isentropic  process,  takes 
;al  gas  correlations  to  determine  the 


Nomenclature 


A 

H 

P 

S 

T 

V 

w - 

I - 

Subscripts 

e - 

f 

g 

o - 

r - 
sat  - 


nozzle  exit  cross  sectional  area 

enthalpy 

pressure 

entropy 

temperature 

specific  volume 

mass  flow  rate 

quality 


exit 

saturated  liquid 
saturated  vapor 
stagnation  (upstream) 
reference 

saturation  temperature 


Mathematical  Description 


The  input  pressures  P^  and  Pg  are  first  converted  from  gage  pressures  to 
absolute  pressures  by  adding  the  second  parameter,  Pj.,  which  represents  the 
absolute  pressure  corresponding  to  zero  gage  pressure.  The  entropy,  inlet 
enthalpy,  saturation  temperature,  saturation  entropy  and  outlet  temperature 
are  then  determined  by  calling  the  steam  property  function  subprograms  de- 
scribed in  chapter  3. 

S = S„  = = SS  ( P„,T„  ) 

=0  = HS  ( P„,T„  ) 


' TSATS  ( P^  ) 


Sg  = SSATS  ( ) 

T = TPSS  < P.  ,S  ) 
e e 

Condensation  occurs  if  Sg  is  greater  than  Sg.  If  condensation  does  occur, 
then  the  concept  of  quality  is  used  to  determine  the  exit  properties. 


Sf  = SSATW  ( Tgg^  ) 
Vf  = VSATW  ( Tgg^  ) 
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V,  “ VSATS  ( ) 

Hf  - HSATW  ( T„t  ) 

Bg  - HSATS  ( ) 

S,  - Sf 

V,  = Vf  + x(  Vj  - Vf  > 

Hg  “ + x(  Hg  “ H£  ) 

If  condensation  does  not  occur,  then  the  outlet  properties  are  determined  as 
follows  : 

V,  = vs  ( P^.T,  ) 

= HS  ( P,,T^  ) 

The  steam  mass  flow  rate  can  now  be  calculated.  The  factor  of  1000  in  the  fol- 
lowing equation  converts  the  enthalpy  from  kj/kg  to  J/kg. 

w = [ 2(  Hq-  Hg  )1000 

"e 

TYPE  23  Component  Configuration 


Inputs 

Description 

1 

To 

- upstream  stagnation  temperature  [ C ] 

2 

Po 

- upstream  stagnation  pressure  [ KPa,  gage  ] 

3 

Pe 

- exit  pressure  [ EPa,  gage  ] 

Outputs 

Description 

1 

Te 

- exit  temperature  [ C ] 

2 

w 

- steam  mass  flow  rate  [ kg/s  ] 

Parameters 

1 

A 

- nozzle  exit  cross  sectional  area  [ m^  ] 

2 

Pr 

- absolute  pressure  corresponding  to  zero  gage  pressure 
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4.24  TYPE  2 4 : IDEAL  GAS  NOZZLE 


General  Description 

The  Type  24  subroutine  models  a convergent  nozzle  for  any  ideal  gas  and 
assumes  the  nozzle  to  be  isentropic,  adiabatic  and  without  condensation.  The 
model  calculates  a mass  flow  rate  and  exit  temperature  given  the  upstream 
stagnation  pressure  and  temperature  and  the  downstream  pressure.  The  model 
uses  absolute  pressures  internally,  but  allows  the  input  and  output  pressures 
to  be  gage  pressures.  The  second  parameter  is  the  absolute  pressure  corre- 
sponding to  zero  gage  pressure.  The  third  and  fourth  parameters  specify 
properties  of  the  ideal  gas. 

Nomenclature 


A 

P 

T 

w 

Y 


nozzle  exit  cross  sectional  area 

pressure 

temperature 

mass  flow  rate 

specific  heat  ratio,  Cp/Cv 


Subscripts 


c - critical 

e - exit 

o - stagnation 


Mathemat ical  Description 


The  two  input  pressures  are  first  converted  from  gage  to  absolute  pressures  by 
addition  of  the  second  parameter,  which  represents  the  absolute  pressure  cor- 
responding to  zero  gage  pressure.  The  critical  pressure  ratio  is  then  defined 
as  : 


Pc  = Cl 


-C2 


where  Cl  = — ^ — 
1 + y 


and  C2  = 


1 - Y 


The  method  used  to  determine  the  flow  rate  depends  on  the  critical  pressure 
ratio  and  on  another  pressure  ratio,  Ra: 

Pe 

Ra  = — 

Po 

The  maximum  flow  rate  is  obtained  when  the  pressure  ratio,  Ra,  is  less  than  or 
equal  to  the  critical  pressure  ratio: 


C2-C3 


w = 1000  A Po  [ ] 
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where  C3 

Y 

Otherwise,  the  flow  rate  is  defined  by 


w = 1000 


ri  *r3  r3  • ^ 

A Po  [ i ] 


Lastly,  the  exit  temperature  is  defined  by  the  following  equation: 


TYPE  24  Component  Configuration 
Inputs  Description 

1 Pq  - upstream  stagnation  pressure  [ KPa,  gage  ] 

2 - downstream  pressure  [ KPa,  gage  ] 

3 ^o  ~ stagnation  temperature  [ C ] 

Outputs  Description 

1 w - mass  flow  rate  [ kg/ s ] 

2 T„  “ exit  temperature  [ C ] 

Parameters  Description 

1 A - nozzle  exit  cross  sectional  area  [ m2  ] 

2 P - absolute  pressure  corresponding  to  zero 

gage  pressure  [ KPa  ] 

3 Y “ ratio  of  specific  heats  [ - ] 

4 R - gas  constant  per  unit  mass  [ J/ (kg-K)  ] 
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4.2  5 TYPE  2 5:  STEAM  TO  AIR  HEAT  EXCHANGER 


General  Description 


This  component  represents  a cross  flow  heat  exchanger  consisting  of  hori- 
zontal tubes  with  circular  fins.  Air  flowing  across  the  finned  tubes  is  heat- 
ed. while  steam  inside  the  tubes  is  condensed.  The  steam  side  is  assumed  to 
be  equipped  with  a trap  which  allows  condensate  but  not  vapor  to  leave  the 
tubes.  No  attempt  is  made  to  compute  the  condensate  temperature  at  the  exit. 
Instead,  the  condensate  outlet  temperature  is  assumed  to  be  a fixed  number  of 
degrees  below  the  saturation  temperature.  In  addition,  the  average  tempera- 
ture of  steam  and  condensate  is  assumed  equal  to  the  saturation  temperature, 
and  condensation  is  assumed  to  occur  at  constant  pressure.  The  method  used  to 
represent  air  temperature  dynamics  assumes  that  the  dynamics  are  due  to 
changes  in  the  steam  pressure  and  thus  in  the  saturation  temperature,  rather 
than  to  changes  in  the  air  inlet  temperature  or  flow  rate. 


Inputs  to  the  model  are  the  steam  pressure  and  inlet  temperature,  the  air 
inlet  temperature  and  flow  rate,  and  the  air  outlet  temperature,  which  is  ob- 
tained from  the  first  output  of  the  model.  Other  outputs  are  the  condensate 
outlet  temperature,  the  mass  flow  rate  of  the  steam  (which  is  determined  by 
the  condensation  rate),  and  the  rate  of  heat  transfer  to  the  air. 


Nomenclature 


A 

a 

Ai 

Cm 

Cpa 

Cpw 

Df 

»i 

Do 

Fm 

g 

l^i 


h'® 

“si 

®wo 

Ka 

^f 

Kt 

Kw 

^tu 

Nu 

Pr 

Ps 

Q 

Q 

ss 

Re 

Rr 


minimum  air  flow  area  (nr) 

tube  inside  surface  area  (m^) 

primary  (tube)  outside  surface  area  (m"^) 

secondary  (fin)  surface  area  (m'^) 

thermal  capacitance  of  heat  exchanger  (kJ/C) 

specific  heat  of  air  [kj/ (kg  C) ] 

specific  heat  of  liquid  water  [kJ/ (kg  C)] 

fin  diameter  (m) 

tube  inside  diameter  (m) 

tube  outside  diameter  (m) 

number  of  fins  per  meter 

acceleration  of  gravity  (m/s^) 

inside  surface  heat  transfer  coefficient  [kW/ (m^  C)] 

outside  surface  heat  transfer  coefficient  [kW/ (m^  C) ] 

latent  heat  of  vaporization  of  water  (kJ/kg) 

steam  enthalpy  at  inlet  (kJ/kg) 

condensate  enthalpy  at  outlet  (kJ/kg) 

thermal  conductivity  of  air  [kW/ (m  C) ] 

fin  thermal  conductivity  [kW/(m  C) ] 

tube  thermal  conductivity  [kW/ (m  C) ] 

thermal  conductivity  of  liquid  water  [kW/(m  C) ] 

number  of  transfer  units  ( - ) 

Nusselt  number,  based  on  tube  outside  diameter  ( - ) 

Prandtl  number  ( - ) 

steam  pressure  (kPa,  absolute) 

dynamic  rate  of  heat  transfer  to  air  (kW) 

steady  state  rate  of  heat  transfer  (kW) 

Reynolds  n\imber,  based  on  tube  outside  diameter  ( - ) 

ratio  of  finned  tube  surface  area  to  unfinned  tube  surface  area 
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- inside  surface  heat  transfer  resistance  (C/kW) 

Rq  - outside  surface  heat  transfer  resistance  (C/kW) 

- ratio  of  inside  surface  heat  transfer  resistance  to  total  heat 
transfer  resistance  ( - ) 

R^  - tube  heat  transfer  resistance  (C/kW) 

Tgj  - air  inlet  temperature  (C) 

- air  outlet  temperature  (C) 

Taos  “ steady  state  air  outlet  temperature  (C) 

Tgg  “ difference  between  saturation  temperature  and  condensate  outlet 
temperature  (C) 

Tgg  - saturation  temperature  of  steam  (C) 

- steam  inlet  temperature  (C) 

“ tube  inside  surface  temperature  (C) 

T^o  ~ condensed  water  outlet  temperature  (C) 
w-  - air  mass  flow  rate  (kg/s) 

Wg  - steam  mass  flow  rate  (kg/s) 

6 - fin  thickness  (m) 

r\  - fin  efficiency  (-) 

- air  viscosity  [kg/ (m  s)] 

- liquid  water  viscosity  [kg/ (m  s)] 

- density  of  water  (kg/m^) 

T - heat  exchanger  time  Constant  (s) 


Mathematical  Description 

The  model  first  determines  the  inlet  steam  enthalpy,  the  saturation 
temperature,  and  the  outlet  temperature  and  enthalpy,  using  property  functions 
described  in  section  3.5: 

Hsi  ” BS(Ps.  Tsj) 


T 

ss 

= TSATS(Pg) 

T 

= T 

- T„.. 

•^wo 

ss 

■^sc 

^wo 

= HSATW(Tss 

Next  the  outside  heat  transfer  coefficient  is  computed,  using  a correlation 
from  Hausen  [35]  for  air  flow  over  finned  tubes: 


Nu  = C ReO.625  g^-0.375  pfl/3 

In  this  equation,  the  constant  C is  0.30  for  in-line  tube  rows  or  0.45  for 
staggered  tube  rows.  A value  of  0.30  is  used  in  the  present  model.  The 
Prandtl  number,  Pr,  is  evaluated  at  the  bulk  air  temperature,  using  air  prop- 
erty functions  described  in  section  3.6.  Rf,  the  ratio  of  outside  surface 
area  to  the  area  of  unfinned  tubes  with  the  same  outside  diameter,  is  calcu- 
lated from  the  parameters.  The  characteristic  length  used  in  the  Reynolds 
number.  Re,  and  in  the  Nusselt  number,  Nu,  is  the  tube  outside  diameter. 
Thus , 


Re  = 

Ma^a 
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and 


ho  = Nu  Ka 
Do 

The  fin  efficiency,  ti , is  determined  by  the  use  of  subroutine  SUFED  described 
in  section  3.3.  The  outside  heat  transfer  resistance  is  then  given  by 

•^o  = t + Ap)  ]"^ 


The  heat  transfer  resistance  of  the  heat  exchanger  tubes  is 

(Dp  - Dj> 

^t  = (2KtAi> 

The  inside  surface  temperature  of  the  pipe,  T^,  and  the  inside  heat  transfer 
coefficient,  h^,  are  determined  iteratively  from  the  following  equations: 

r F tO.25 

hj  = 0.612  [ \ Pw  ®5lg.  _ 1 

MwDidss-  Tt)  ■' 


^ss  ^r^^ss  ^a^ 

Rj 

Ri  + Rt  + Rq 

Ri  = (h-A^)"^ 
and  the  correlation  for  h^  is  from  Kern  [36]. 

The  steady  state  rate  of  heat  transfer  is  used  to  find  the  steady  state  air 
outlet  temperature  and  the  rate  of  condensation,  which  is  equal  to  the  enter- 
ing steam  mass  flow  rate: 


T.  = 


where 


R_  = 


®ss  = kiAid,,  - Tt) 


T = T • + 
■^aos  ai 


Q 


s s 


'’'a'pa* 


Q 


s s 


^Hsi  - H^o) 


Air  temperature  dynamics  are  computed  by  the  approximate  method  of  Myers  et 
al.  [37],  assuming  that  the  thermal  capacitance  of  the  heat  exchanger  walls  is 
much  greater  than  that  of  the  air. 

dT  T - T 
= aos ao 

dt  T 
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where 


C« 

«».Cp, 


2 ezp(-z)  tinh(z) 


R^{[sinh(z)  + cosh( z)  ]ezp(*-z)  - exp(-N^^)) 


N 


z = 


tu 


2(1  - Ry) 


-1 


Finally,  the  dynamic  rate  of  heat  transfer  to  the  air  is  calculated: 
Q = WgCpg^(Tmj  ~ 

TYPE  25  Component  Configuration 


Inputs 

1 

2 

3 

4 

5 

Outputs 

1 

2 

3 

4 


Description 

Pg  - steam  pressure  [ KPa,  absolute  ] 

Tg£  “ inlet  steam  temperature  [ C ] 

T^£  - inlet  air  temperature  [ C ] 

~ outlet  air  temperature  [ C ) 

- mass  flow  rate  of  air  [ kg/s  ] 
Description 

- outlet  air  temperature 

- outlet  temperature  of  condensate  [ C ] 

- mass  flow  rate  of  steam  [ kg/ s ] 

Q - rate  of  heat  transfer  to  air  [ kW  ] 


a 

T 

ao 

^s 
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Parameters 


Description 


1 

- subcooling  of  condensate  below  saturation  [ C ] 

2 A, 

- minimum  air  flow  area  [ m2  ] 

3 Ai 

- inside  heat  transfer  area  [ m^  ] 

4 Ap 

A 

- primary  (tube)  outside  surface  area  [ m“^  ] 

5 A^ 

- secondary  (fin)  surface  area  [ m*^  ] 

6 Di 

- tube  inside  diameter  [ m ] 

7 *>o 

- tube  outside  diameter  [ m ] 

8 K^ 

- tube  thermal  conductivity  [ kJ/ (kg  K)  ] 

9 Dj 

- fin  diameter  [ m ] 

10  Kf 

- fin  thermal  conductivity  [ kJ/ (kg  K)  ) 

11  6 

- fin  thickness  [ m ] 

12  Fm 

- number  of  fins  per  meter 

13 

ni 

- thermal  capacitance  of  metal  (fins  and  tubes)  [ kJ/K] 
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4.26  TYPE  26:  CONTROL  SIGNAL  INVERTER 


General  Description 

This  component  is  a very  simple  representation  of  an  inverting  relay. 
Its  single  input  is  a control  signal  between  zero  and  one,  and  its  output  is 
one  minus  the  input. 

In  a typical  simulation,  a single  control  variable  might  represent  both 
the  output  of  a controller  and  the  relative  position  of  a valve  or  damper. 
The  TYPE  8 and  TYPE  16  controller  models  operate  on  the  assumption  that 
increasing  the  output  control  signal  will  increase  the  controlled  variable. 
Depending  on  the  system,  it  may  be  necessary  to  use  a TYPE  26  inverter  between 
a controller  and  a valve  or  damper  in  order  to  maintain  the  validity  of  this 
assumption. 

Nomenclature 

C^  - input  control  signal 

Cq  - output  control  signal 

Mathematical  Description 


0*  < Co  < 1. 

TYPE  26  Component  Configuration 


Input 


Description 


1 


- input  control  signal 


Output 


Description 


1 


Co  “ output  control  signal 


Parameters:  NONE 
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